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We present a new three-dimensional general-relativistic hydrodynamic evolution scheme coupled 
to dynamical spacetime evolutions which is capable of efficiently simulating stellar collapse, isolated 
neutron stars, black hole formation, and binary neutron star coalescence. We make use of a set 
of adapted curvi-linear grids (multipatches) coupled with flux-conservative cell-centered adaptive 
mesh refinement. This allows us to significantly enlarge our computational domains while still 
maintaining high resolution in the gravitational-wave extraction zone, the exterior layers of a star, or 
the region of mass ejection in merging neutron stars. The fluid is evolved with a high-resolution shock 
capturing finite volume scheme, while the spacetime geometry is evolved using fourth-order finite 
differences. We employ a multi-rate Runge-Kutta time integration scheme for efficiency, evolving 
the fluid with second-order and the spacetime geometry with fourth-order integration, respectively. 
We validate our code by a number of benchmark problems: a rotating stellar collapse model, an 
excited neutron star, neutron star collapse to a black hole, and binary neutron star coalescence. The 
test problems, especially the latter, greatly beneflt from higher resolution in the gravitational-wave 
extraction zone, causally disconnected outer boundaries, and application of Cauchy-characteristic 
gravitational-wave extraction. We show that we are able to extract convergent gravitational-wave 
modes up to {£,m) = (6,6). This study paves the way for more realistic and detailed studies of 
compact objects and stellar collapse in full three dimensions and in large computational domains. 
The multipatch infrastructure and the improvements to mesh refinement and hydrodynamics codes 
discussed in this paper will be made available as part of the open-source Einstein Toolkit. 

PACS numbers: 04.25.D-, 04.30.Db, 97.60. Bw, 02.70.Bf, 02.70.Hm 



I. INTRODUCTION 

Some of the most interesting relativistic astrophysi- 
cal phenomena such as stellar collapse, black hole forma- 
tion, or binary neutron star coalescence, require numeri- 
cal simulations on large computational domains, involve 
many different length scales, and are intrinsically three- 
dimensional (3D). Due to their extreme nature in terms 
of fluid densities and velocities, an accurate treatment 
of general-relativistic (GR) gravity is required. Depend- 
ing on the problem, magnetic field evolution and neu- 
trino interactions may also be required. Thus, numerical 
computations in relativistic astrophysics are truly multi- 
physics, and as such, are especially demanding in terms 
of computational modeling technology and resources. 

Current state of the art 3D GR hydrodynamic simula- 
tions in the context of stellar collapse [lj-i5j or binary neu- 
tron star coalescence [BHIB] (see [2] for a recent review) 
are based on Cartesian grids with adaptive mesh refine- 
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ment (AMR) . As the domain is enlarged or the resolution 
increased, such grids pose a serious bottleneck in terms 
of the computational power that is required, even with 
AMR. Since Cartesian grids scale as in terms of the 
number N of grid points along one spatial direction in 3D, 
available computational resources are rapidly exhausted 
when additional points in each coordinate direction are 
added. The symmetry of the computational problem, 
however, is essentially spherical, at least at some distance 
from the central region of the simulation. Thus, Carte- 
sian grids are wasteful with respect to angular resolution 
when the problem becomes symmetrically spherical. 

For instance, stellar collapse proceeds in approximately 
spherical or axisymmetric terms (e.g. |15tfl8] ). At later 
times, various hydrodynamic instabilities (e.g. convection 
and instabilities of the shock) break this symmetry. The 
global features, however, remain approximately spherical 
or axisymmetric. 

In the case of coalescing binary neutron stars, the 
central region containing the two neutron stars is not 
of spherical symmetry. At larger distances and in the 
gravitational- wave (GW) zone, however, the problem 
becomes spherical. The gravitational-wave extraction 
zone must generally be located at large radii in order 
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to limit near-zone effects in the extracted wave. But 
even with more sophisticated techniques such as Cauchy- 
characteristic extraction |19ti2^ that allow us to ex- 
tract gauge-invariant GWs at future null infinity ^7+, it 
is necessary to enlarge the domain sufficiently so that 
constraint- violating modes generated at the outer bound- 
ary are causally disconnected from the interior evolution 
and the wave-extraction zone. These constraint- violating 
modes are generated due to the lack of constraint- 
preserving outer boundary conditions for the Einstein 
equations (see [53] for a recent review) for certain types 
of evolution systems (including the common BSSN sys- 
tem), and travel at the speed of light P^[?7] to the inte- 
rior of the domain. Without these systematic errors, the 
evolution and wave extraction would generally be more 
accurate. Furthermore, in case mass is ejected during 
and after merger, enlarging refinement levels to track the 
evolution of the ejected material becomes very expensive. 

It therefore seems natural to apply spherical grids to 
maintain high resolution also in the outer regions of the 
domain. The computational effort when using spherical 
grids scales linearly with the number of radial points N, 
assuming constant angular resolution. Thus, spherical 
grids can give a tremendous performance improvement 
when the domain is enlarged or the (radial) resolution 
increased. 

Spherical grids have been widely used for many astro- 
physical problems, including stellar collapse (e.g., |281 - 
l32]). core-collapse supernovae (e.g., |33H35] ). oscillations 
of neutron stars (e.g., [Ml 132]), neutron star magneto- 
sphere (e.g., [33), accretion onto black holes [31], and 
simulations of accretion disks (e.g., |401 - H5] ) . Unfortu- 
nately, the standard spherical-polar coordinate system 
imposes a serious difficulty along the axis and the poles, 
where special care must be taken to regularize the fields 
and to numerically evolve them |43fH5] . But even with 
proper regularization applied, the angular and radial dis- 
tribution of grid points is non-optimal in the sense that 
they cluster at the poles and at the coordinate origin. In 
addition, spherical grids are less suited in regions where 
the underlying symmetry is non-spherical, e.g., in the 
vicinity of a binary neutron star system, or the highly tur- 
bulent and convective region behind the accretion shock 
in a core-collapse supernova. 

In order to handle multiple regions of different symme- 
try within the same simulation, multipatch (sometimes 
also called multiblock) schemes have been developed for 
a wide range of physics and engineering applications. 
The idea is to cover the simulation domain with mul- 
tiple curvi-linear coordinate "patches". Each patch is lo- 
cally uniform. Diffeomorphic mappings from local to the 
global coordinates enable to represent a wide range of 
grid shapes in different regions of the simulation. One 
such example is given in Fig. [T] In this setup, a cen- 
tral Cartesian patch is surrounded by six "inflated cube" 
spherical grid patches. This is a natural configuration for 
our purposes. The aspherical region of a collapsing star 
or a merging binary is best modeled by a central Carte- 



sian patch, capable of AMR. The gravitational-wave zone 
and/or the outer layers of a star are best modeled by 
the more efficient spherical grids. This allows us to em- 
ploy large domains at high resolution with modest com- 
putational cost. Notably, the outer boundary can be 
causally disconnected from the interior evolution and the 
gravitational- wave extraction zone. 

Within the context of numerical relativity and rela- 
tivistic astrophysics, multipatch schemes have already 
been successfully applied in a range of different problems 
ranging from simulations of accretion disks |47l I48j , hori- 
zon finding [39], wave extraction [SD], single black holes 
[511 [S2] , orbiting black holes [33^ , relativistic fluid evolu- 
tions on fixed backgrounds [54 , elliptic and initial data 
solvers |55fl55] . to characteristic evolutions of Einstein's 
equations |20l I5U1 [M] . Multidomain spectral methods 
have been successfully applied to vacuum binary black 
hole evolutions yielding high accuracy and efficiency [62]- 
I67| using a dual-coordinate frame method [68J . The same 
multidomain spectral code SpEC, coupled to a finite vol- 
ume fluid solver, has also been used to simulate neutron 
star black hole mergers |69fl7^ . Neither of the works 
above, however, make use of AMR for the fluid flelds, and 
thus are limited in the respective range of astrophysical 
applications. In particular, efficient simulations of stellar 
collapse and black hole formation require AMR in the 
central region of the collapsing star. Also, the near-fleld 
region in simulations of binary neutron star coalescence 
substantially beneflt from AMR, in particular when ma- 
terial is ejected in the post-merger phase. 

In the context of vacuum binary black hole merger sim- 
ulations, multipatch schemes combined with AMR have 
been successfully applied |73ff79] . We base our code on 
the Llcona infrastructure developed in |73j . which makes 
use of the Cactus computational toolkit [50] and the 
Carpet AMR driver [HJIHll. We extend the original pure 
vacuum scheme to include full matter dynamics using 
the publicly available GR hydrodynamics code GRHydro, 
which is part of the EinsteinToolkit [83,. We thus 
present the flrst successful multipatch scheme capable of 
AMR that can stably evolve fluid dynamics coupled to 
fully GR spacetime dynamics. 

In addition, we make a number of improvements: (i) 
We extend the AMR driver Carpet to support cell- 
centered mesh reflnement, which allows us to apply re- 
fluxing, a technique to maintain conservation of mass, 
energy and momentum fluxes across mesh refinement 
boundaries [84] (see [85] for a recent application to GR 
hydrodynamics). This greatly improves conservation of 
mass in our simulations of stellar collapse, especially in 
the postbounce evolution, (ii) We apply enhanced PPM 
(piecewise parabolic method) reconstruction |86) [87] . 
which significantly improves the numerical accuracy and 
the behavior of the constraints, (iii) To improve the 
execution speed of the simulations, we apply multirate 
Runge-Kutta (RK) time integration (e.g. [88] l89] ) in 
which the spacetime is evolved with a standard fourth- 
order RK method, whereas the fluid is evolved with a 
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second order RK scheme without significant loss of accu- 
racy. This reduces the number of intermediate steps in 
the fluid evolution, which dominates in terms of proces- 
sor cycles compared to spacetime evolution, in particular 
when using a microphysical equation of state. 

We apply the new code to a number of benchmark 
problems, including the evolution of a single isolated and 
perturbed neutron star, the collapse of a rotating stel- 
lar core, the collapse of a neutron star to a black hole, 
and the merger of a binary neutron star system. We 
investigate the accuracy and convergence of each test 
problem. This is an important code verification towards 
our program to carry out fully 3D simulations of core- 
collapse supernovae (see [90 for a recent application of 
our scheme) and black hole formation in the context of 
the coUapsar scenario for long gamma-ray bursts. The 
new multipatch scheme allows us to significantly enlarge 
the computational domain by maintaining a fixed angular 
resolution. This is useful in many ways: (i) we are able 
to causally disconnect the outer boundary from the inte- 
rior evolution and the gravitational-wave extraction zone, 
thus avoiding systematic errors from the approximate 
and non-constraint preserving artificial outer boundary 
condition, (ii) we have a larger wave-extraction zone with 
higher overall resolution, thus making it possible to ex- 
tract higher-order than the dominant GW modes, (iii) 
in binary neutron star mergers, ejected material can be 
tracked out to large radii with relatively high resolution, 
(iv) the number of mesh refinement levels can be de- 
creased, leading to better parallel scaling. As a result, 
our multipatch scheme can efficiently evolve models of 
stellar collapse in full 3D (see also fSUJ), and is capable 
of more accurate gravitational-wave extraction in models 
of binary neutron star mergers. In the latter test prob- 
lem, we extract convergent gravitational-wave modes up 
to (£,m) = (6,6). 

This paper is organized as follows. In Sec. |II A| and 
|IIB[ we first review the underlying hydrodynamic and 
spacetime evolution systems and how we solve them nu- 
merically. Subsequently, in Sec. II C we present our ap- 



proach to multipatches and their numerical implementa- 
tion. We also discuss our implementation of cell-centered 



AMR (Sec. IIDl, and describe multirate RK time inte- 
gration (Sec. HE). Finally, in Sec. Ill we present de- 



tailed tests of isolated perturbed and unperturbed neu- 
tron stars, collapsing stellar cores, neutron star collapse 
to a black hole, and merging binary neutron stars. We 
conclude and summarize our findings in Sec. |IV[ In an 
appendix, we presents basic tests with shock tubes (Ap- 
pendix |A|, we review the enhanced PPM scheme as de- 
veloped in |86[ 157] (Appendix |B]) , discuss our treatment 
of the artificial low-density atmosphere (Appendix [C| , 
present an optimized ghost-zone update scheme to im- 
prove the parallel scaling (Appendix |d|, describe our 
volume integration scheme for overlapping grids (Ap- 
pendix |e]) , and investigate the infiuence of boundary ef- 
fects on binary neutron star merger dynamics and wave 
extraction (Appendix [f| . 



II. METHODS 

A. General-Relativistic Hydrodynamics 

We base our code on the open-source GR hydrodynam- 
ics code GRHydro that is part of the EinsteinToolkit 
inj and is described in pT] [551 [52] . 

We introduce primitive variables in the form of the 
fiuid density p, the fiuid's specific internal energy e, and 
the fiuid 3-velocity as seen by Eulerian observers at rest 
in the current spatial 3-hypersurface | i93j . 
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where is the fiuid 4- velocity, W = {1 — v^Vi)' 
the Lorentz factor, and a and /3* are lapse and shift. 



respectively (to be introduced in Sec. II B I. In terms of 



the 3-velocity, the contravariant 4-velocity is then given 

by 
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and the covariant 4-velocity is 

uo = Wiv'Pi - a) , u, = Wvi 



(2) 



(3) 



The evolution equations are written in the Valencia 
form of GR hydrodynamics [Ml [SS] as a first-order hy- 
perbolic fiux-conservative evolution system for the con- 
served variables D, 5*, and r which are defined in terms 
of the primitive variables p, e, w*. 



D - ^pW, 
S' = ^phW^v\ 
T = ^{phW^-P)-D, 



(4) 



where 7 is the determinant of the 3-metric (see 
Sec. II B I, and the quantities P, and h = 1 + e -\- P/p 



denote pressure, and specific enthalpy, respectively. The 
evolution system then becomes 
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with 
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S = 



a 



9 In a 

dxf^ 



(6) 



Here, — — P^/a, F^^, are the 4-Christoffel sym- 
bols, and T^^^ is the stress-energy tensor. The pres- 
sure P = P(p, e, {X;}) is obtained via our equation of 
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state module, which is capable of handling a set of dif- 
ferent equations of state, including microphysical finite- 
temperature variants. The {Xi} are additional composi- 
tional variables of the matter such as the electron fraction 
Ye, which are used for microphysical equations of state. 
In the present work, however, we resort to simple (piece- 
wise) polytropic and ideal gas (F-law) equations of state. 

The above evolution equations are spatially discretized 
by means of a high-resolution shock-capturing (HRSC) 
scheme using a second-order accurate finite- volume algo- 
rithm. The equations are kept in semi-discrete form and 
first-order (in space) Riemann problems are solved at cell 
interfaces with the approximate HLLE solver |Mj . 

The states at cell interfaces are reconstructed using 
a new and improved variant of the piecewise parabolic 
method (PPM) ^ [S71 [97]. As noted in [Ml [HZ], the 
original PPM scheme [HZ] has the side-effect of flatten- 
ing local smooth extrema which are physical, thus lim- 
iting the accuracy. In the present context of simulating 
compact objects, one naturally has extrema at the stel- 
lar center(s) where the matter density is largest. We 
find that the original PPM scheme reduces the accuracy 
there, which then strongly affects the overall accuracy of 
our simulations (see Sec. Ill and also Fig. 251. Ref. [86j, 
further refined by Ref. [87j, suggests modifications to the 
original limiter which can distinguish between smooth 
maxima that are part of the solution, and artificial max- 
ima that may be introduced at shocks and other discon- 
tinuities. While smooth maxima need to be retained as 
part of the solution, artificial maxima must be avoided to 
suppress Gibbs phenomenon at shocks and other discon- 
tinuities. We summarize the procedure for "enhanced" 
PPM reconstruction in Appendix [B| 

We note that under certain conditions, the require- 
ment that the modulus of the reconstructed primitive 
velocity must stay below the speed of light c may be 
violated. This can happen, since the primitive veloc- 
ity is a bounded function (bounded by the requirement 
Viv'^ < c^), and the enhanced PPM reconstruction scheme 
does not enforce this constraint close to any occuring 
extrema. Thus, the enhanced PPM scheme may recon- 
struct velocity components that result in a velocity mod- 
ulus equal to or slightly larger than the speed of light 
near extrema. To avoid this problem, we reconstruct 
Wv^, i.e. the Lorentz factor W times the primitive ve- 
locity . The quantity Wv^ is unbounded and thus does 
not require special treatment near extrema. 

The time integration and coupling with curvature 
(Sec. 



II B I are carried out with the Method of Lines 



(see Sec7[IIEl. 

After each evolution step, we compute the primitive 
quantities from the evolved conserved quantities. Since 
the primitive quantities are implicit functions of the con- 
served ones, it is necessary to use a numerical root finding 
algorithm. As described in, e.g. |53], this is done via a 
Newton-Raphson scheme. 

In some rare situations, the initial guesses for the 
root finding procedure are not well-posed, and cause the 



Newton-Raphson scheme to fail to converge. In partic- 
ular, we find this behavior at the surface of a neutron 
star, when the latter is threaded by an AMR boundary 
and refluxing is active. In this case, we resort to a simple 
bisection algorithm which converges more slowly, but is 
more robust. 

In regions of the computational domain, where we have 
physical vacuum, we employ an artificial low density "at- 
mosphere" (see Appendix [C]) . In order to reduce the in- 
fluence of the artificial atmosphere on the curvature evo- 
lution, we exponentially damp the stress-energy tensor 
Tf^i, to zero outside a given radius. More specifically, 
we introduce the radius dependent stress-energy damp- 
ing T^^ — >■ X{r)T^i, with the damping factor 



Mr) = <! l(l-tanh(«r^i-^) 




for r < i?o, 
otherwise, (7) 
for r > Ri, 



where the damping is applied between the two radii Rq < 
R,. 

At outer boundaries, we apply a copy-from-neighbor 
(flat) boundary condition for the evolved fluid quantities. 

Finally, in order to be compatible with multipatch dis- 
cretization, we need to introduce additional coordinate 
transformations as described in Sec. Ill C~3l below. 



B. Curvature Evolution 

The spacetime evolution is performed by a variant of 
the BSSN evolution system |5^102| and is implemented 
in the CTGcumna curvature evolution code [73J , which was 
developed for arbitrary coordinate systems mapping the 
spatial domain. 

The standard BSSN system is derived from a 3-f 1 split 
of spacetime resulting in a foliation in terms of spatial 
hypersurfaces along a timelike vector field. It introduces 
the following set of evolved variables 



<j), 7afc, K, Aab, F", 

which are solved according to 

dt^^~laK+ld,(3\ 
6 6 



(8) 



(9a) 



dtjab = - 2aAab + 13'daab + 27.(a5b)/3' (9b) 

dtK = - D.D'a + a{AijA'^ + \k^) + (5'd,K (9c) 

+ Aira (pADM + S) , 
dtAab ^e-^^i-DaDba + aR^b)'^^ + p^d^A^b (9d) 

+ 2i,(,9b)/3' - "^Aabd^P' 



iire-^'t'a (S, 



ab) 



5 



dtt- + lr'^^^J^' - ra,/3" (9e) 

3 

- 167ra7"S'„ 

where Da is the covariant derivative determined by 
the conformal 3-metric 7^6, and "TF" indicates that the 
trace-free part of the bracketed term is used. 

Above, we show the "</>"- variant of the BSSN system. 
Our curvature evolution code also provides the "x"- and 
" H^"- variants of the evolution system (see [73^ for details) . 
Here, we employ the 0- variant. 

The stress-energy tensor T^^ is incorporated via the 
projections 

PADM ^ (Too - 2/3*ro, + (3^/3^^^) , (10) 
S := rr,,, (11) 

Sa ■■= -^{Toa-f3'Ta,) , (12) 

{Sabf^ ■■= (Tab - ^e4^S7ab) • (13) 

After each evolution step, the evolved curvature vari- 
ables ([8| are transformed (via an algebraic relation) 
to the standard ADM variables {gij,Kij} (e.g., |103| 1. 
where gij is the (physical) 3-metric, and Kij the extrin- 
sic curvature. The ADM variables are used to couple 
the curvature evolution to the hydrodynamic evolution 
scheme, i.e., our hydrodynamic scheme uses the physical 
3-metric gij rather than the evolved conformal 3-metric 
7ab above. 

The lapse gauge scalar a is evolved using the 1 -I- log 
condition |104| . 

dta - f3'd^a = -2aK, (14) 

while the shift gauge vector /3° is evolved using the hy- 
perbolic F-driver equation |105j . 

dtl3'' - p'd,^'' = ^B" , (15a) 

dtB" - l3^djB' = 9tf ° - /3'a,f ° - q{r)riB'' , (15b) 

where 77 is a parameter which acts as a (mass dependent) 
damping coefficient. To avoid certain stability issues with 
the gauge arising in the far-field regime [106 , the damp- 
ing coefficient is allowed to spatially change, either by 
some dynamic evolution [107', or by a fixed prescription. 
We use the simple prescription for a radial fall-off of i] 
given in |106| If not stated otherwise, we use a fall-off 
radius of i? = 250 M©. 

The 3+1 decomposition of the Einstein equations also 
results in a set of constraint equations. The Hamiltonian 
constraint equation reads 

H = + ^2 - K^jK'^ - 167rpADM = , (16) 



where R^^^ denotes the 3-Ricci scalar, and the momentum 
constraint equations read 

M" EE ~ 7"^) ~ SttS"" = . (17) 

We do not actively enforce the constraints during evolu- 
tion, but rather check how well our numerically obtained 
metric quantities satisfy the constraints over the course 
of the evolution. Thus, this offers a valuable accuracy 
monitor for the curvature evolution. 

The spacetime equations are discretized using fourth- 
order finite difference operators |108j . The finite differ- 
ence stencils are centered. An exception are the advec- 
tion terms of the form /3'9i, which use operators that are 
upwinded by one stencil point towards the local direction 
of the shift vector [73J. 

Consistent with the order of accuracy of spatial finite 
difference derivatives, we also apply Kreiss-Oliger dissi- 
pation |108j which is of one order higher than the spatial 
discretization order. In the case of fourth-order differ- 
encing, we thus apply fifth-order dissipation operators. 
Dissipation is added to the right-hand-sides (RHS) of 
the curvature evolution quantities at any time integra- 
tion substep. The strength of the dissipation can be con- 
trolled by a parameter ediss G [0,1]- Unless otherwise 
specified, we use Cdiss = 0.1 throughout this work. 

At outer boundaries, we impose a simple approximate 
radiative boundary condition as described in [73 . Since 
data from this condition are not strictly constraint sat- 
isfying, constraint violating modes are generated at the 
boundary, and travel with the speed of light |261 [27] to 
the interior of the domain where they introduce a sys- 
tematic error in the curvature evolution. 



C. Multipatches 

We build our code on the Llsuna infrastructure de- 
scribed in detail in |73| . This infrastructure implements 
multipatches via an arbitrary number of curvi-linear 
overlapping grid patches using fourth-order Lagrange and 
second-order essentially non-oscillatory (ENO) interpola- 
tion for exchanging data in inter-patch ghost zones be- 
tween neighboring patches. In |72|, only the pure vacuum 
problem was considered. Here, we extend the multipatch 
evolution scheme to include matter. 



1. Patch Systems 

A useful patch system is shown in Fig. ^ the central 
Cartesian patch is surrounded by six spherical infiated- 
cube patches. The nominal^ grids of the spherical 
patches have inner radius i?g, outer radius Rb, radial 



We define the nominal grid as the unique set of points covering 
the entire computational domain, i.e. the nominal grid of a single 



spacing Ai?i, which is allowed to stretch to Ai?2 within 
some finite region, and angular resolution (Ap, Act) per 
angular direction (p, cr). Note that the angles (p, cr) used 
to define the local coordinates of each infiated-cube patch 
do not coincide with standard spherical-polar coordinates 
(see below). The central patch contains a hierarchy of 
refined regions, allowing to place resolution where neces- 
sary. This patch system is particularly useful in problems 
with spherical symmetry at some radius from the central 
source. 

Each grid patch defines local uniform coordinates 
{u, V, w) related to the global Cartesian (x, y, z) coordi- 
nate space by a diffeomorphic relation. For the central 
Cartesian patch depicted in Fig. [T] this relation is triv- 
ially given by the identity function. The infiated-cube co- 
ordinates, however, are defined by non-trivial coordinate 
functions. For each angular patch, we define local an- 
gular coordinates (p, cr) that range over (— 7r/4, +7r/4) x 
(— 7r/4, +tt/A) and can be related to global angular coor- 
dinates (/X, 1/, (j)) (see Fig. [l]) which are given by 

p = rotation angle about the x-axis = arctan(y/z), 

(18a) 

v = rotation angle about the y-axis = arctan(a;/z), 

(18b) 

(j) = rotation angle about the z-axis = aTctan(y / x) . 

(18c) 



rAi?i 



IAR2 




(Ap,Ao-) 



For each angular patch, we have two unique angles (p, a) 
out of the three global angles (p, v, (f>) that parametrize 
the local coordinates. For instance, for the patch normal 
to the positive x-direction, we select 



p = V = arctan(z/2:), 
a = (p = arctan(y/x), 
R = fir), 



(19a) 
(19b) 
(19c) 



where r — y^x'^ + + z^. Similarly, the coordinates 
of the patch along the positive y and z axes are 
parametrized by (p, cr) = (p, </>) and (p, u) = (p, i^), re- 
spectively. The remaining three patches along the nega- 
tive axes are related in a similar way. 

In the radial coordinate direction, we apply radial 
stretching with an appropriate stretching function R = 
f{r). In the stretching region, the physical coordinate 
radius is stretched, corresponding to a smooth decrease 
in radial resolution from spacing Ai?i to spacing Ai?2- 
Outside the stretching region, we keep the radial spacing 
constant. Details can be found in [73) . 



patch excludes ghost points (and additional overlap points; see 
further below) that are shared with a neighboring patch. 



FIG. 1: Depiction of a typical patch system used in our sim- 
ulations. The upper figure schematically shows a. z = slice 
of the employed grids: a central Cartesian grid (patch 0) is 
surrounded by spherical inflated- cube grid patches (patches 
1 — 4 out of a total of six spherical patches). The central grid 
is capable of AMR allowing to refine the resolution at the 
central region of, e.g., a star where the density and curvature 
gradients become large. _Rb and _Rs denote the radii of the 
outer computational boundary and of the boundary between 
spherical and Cartesian grids, respectively. The spherical grid 
has a fixed angular resolution denoted by (Ap, Acr), while the 
radial resolution is allowed to stretch from radial resolution 
ARi to AR2. The lower figure shows a radial R = const, shell 
of the outer spherical grid, comprised of six inflated-cube grid 
patches. Angular points can be uniquely determined by two 
out of three angular coordinates (p, (18 1. Interpolation 
at patch boundaries reduces to ID interpolation. Points are 
almost uniformly distributed across the sphere. 



2. Spacetime Evolution Scheme 

Here, and as described in [73J, the spacetime evolu- 
tion is solved in the global Cartesian (x, y, z) tensor ba- 
sis, where the grid patches are generally distorted, i.e., 
they are not uniform. Derivatives are approximated via 
finite differences in the local coordinate system {u,v,w) 
of each grid patch, where, as required by our finite dif- 
ference scheme, the grid patches are uniform. In order 
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to transform to the global tensor basis, Jacobian trans- 
formations of the form J*j = du^ /dx^ are applied to the 
first and second derivatives at each point, 

_d_^f du£\ _d_^ 

dxi \ dxj J duj ' 

dxidxj \dxidxj J du\ \dxi dxj J dukdui ' 

(20b) 

thus obtaining the derivatives in the global [x, y, z) co- 
ordinate space. The Jacobians are precomputed at each 
grid point. The main advantage of solving the equations 
in the global (x, y, z) basis is simplicity. There is no need 
for inter-patch coordinate basis transformations. Per- 
haps more importantly, the existing code infrastructure, 
and especially analysis tools, do not need to be changed, 
since the assumption of a global Cartesian tensor basis is 
still maintained. 



TABLE I: Required quantities for the hydrodynamic evo- 
lution scheme and their coordinate bases. A tilde denotes 
quantities that need to be obtained by applying a Jacobian 
transformation. The last four quantities are only required for 



microphysical equations 


of state. 






C2uantity 


iype 


global 


local 


metric 


tensor 




9ii 


extrinsic curvature 


tensor 






shift 


vector 






lapse 


scalar 


a 


a 


prim, density 


scalar 


P 


P 


specific internal energy 


scalar 


e 


e 


prim, velocity 


vector 






cons, density 


densitized scalar 




D 


cons, internal energy 


densitized scalar 




T 


momentum 


densitized vector 






stress-energy tensor 


tensor 






Lorentz factor 


scalar 


w 




pressure 


scalar 


p 




prim, electron fraction 


scalar 




Ye 


cons, electron fraction 


densitized scalar 




■(/■con 
^ e 


temperature 


scalar 


T 


T 


entropy 


scalar 


s 


s 



3. Hydrodynamic Evolution Scheme 

Finite volume schemes work well on general unstruc- 
tured meshes. The original implementation of the hy- 
drodynamic evolution code GRHydro, however, assumes 
uniform coordinates. Without a major rewrite of the 
code, we can keep our original scheme by solving the 
Riemann problem in the local frame, where the coordi- 
nates are uniform. This requires no changes to the core of 
the scheme. Any computation simply carries over to the 
local coordinate basis. Effectively, this means that the 
primitive and conserved quantities are thus represented 
in the local coordinate basis. 

Special attention is required when coupling the hydro- 
dynamics solver to the metric solver (Sec. IIC2I. The 



metric solver explicitly computes the metric components 
in the global frame and is thus generally incompatible 
with the hydrodynamic quantities defined in the local 
frame. We therefore introduce the additional step of 
transforming the metric components to the local basis 
before each hydrodynamic RHS step. Correspondingly, 
after each hydrodynamic step, we need to compute the 
stress-energy tensor T^^'^ in the global basis as required 
by the metric solver. 

Since the various analysis tools explicitly assume a 
global coordinate frame for the primitive variables, we 
introduce a separate set of global primitive variables. Ef- 
fectively, this only requires extra memory for the primi- 
tive 3- velocity since the primitive density p and r 
are scalars. Once the primitive quantities are known in 
the global frame, the stress-energy tensor can be directly 
computed in the global frame. 

For clarity, we list the various quantities in their cor- 
responding available coordinate basis in Table [T] 




FIG. 2: Depiction of the second-order ENO inter-patch inter- 
polation scheme used for the fluid variables between two over- 
lapping patches p and q. The inter-patch boundary is indi- 
cated by the vertical line. Each interpolated point in the ghost 
zones (empty boxes) is obtained from an interpolation poly- 
nomial whose stencil is selected based on the local smoothness 
of the interpolated quantity. There are three possible choices: 
left (L) stencil using blue and green points, right {R) stencil 
using green and red points, and first-order (/) stencil using 
only green points. Since none of the stencil points on p are 
allowed to be inter-patch boundary points of p, we need to in- 
troduce a certain number of additional overlap points (filled 
boxes) to ensure that this is the case. 



^. Inter-Patch Interpolation and Coordinate 
Transformation 

Data in the ghost zones of a given grid patch are ex- 
changed via high-order Lagrange polynomial interpola- 
tion for those quantities that are smooth (such as the cur- 
vature evolution variables), and optionally second-order 
essentially non-oscillatory (ENO) interpolation |109| for 
those variables that may contain discontinuities (such as 
the hydrodynamic evolution variables). The scheme is 
depicted in Fig. [2] Ghost points (indicated by empty 
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global 



the "local-to-local" transformation 




FIG. 3: Coordinate systems and their transformations. Local 
coordinates u^pj and of patches p and q, respectively, are 
related via "local-to-local" transformations. "Local-to-local" 
transforms are necessary for fluid variable inter-patch inter- 
polation. The global Cartesian coordinates are used to 
represent the curvature variables and to carry out any analy- 
sis on the curvature or fluid variables, such as gravitational- 
wave extraction, or fluid density oscillation mode analysis. 
Therefore, "global-to-local" and "local-to-global" transforms 
are necessary. 



det 



du 



(p) 



du-', 



(q) 



D 



(p) 



(21) 



between local coordinates M(p) of patch p and local coor- 
dinates (w(q)) of patch q. Hence, after having obtained 
its interpolated value in the "old" basis defined by the lo- 
cal coordinates of patch p, we need to represent it in the 
"new" basis defined by the local coordinates of patch q ac- 



cording to transformation (21 1, before we assign its trans- 
formed value to one of the ghost points of q. Similarly, 
we also need to transform the conserved 3-momentum, 
which transforms as a densitized contravariant vector ac- 
cording to 



det 



^"(P) 



du 



(q) 



du-'. 



(q) 



S}, 



C"(p) 



(22) 



The various coordinate transformations that are re- 
quired in our code are depicted in Fig. [3j 



boxes) on some patch p must be interpolated from points 
from a neighboring overlapping patch q. The inter-patch 
boundary is indicated by a vertical line. For Lagrange in- 
terpolation, in order to maintain maximal accuracy, we 
center the interpolation stencils around the interpolation 
point. The ENO operator, on the other hand, is allowed 
to use second-order off-centered Lagrange interpolation 
stencils according to the local smoothness of the interpo- 
lated fields |109] . In addition, we check if the interpolant 
introduces a local maximum and switch to first order in 
that case. In order to speed up the computation, we pre- 
compute and store all possible stencil configurations for 
each inter-patch ghost point. 

To yield a consistent boundary treatment, we have to 
ensure that an interpolation stencil does not contain any 
ghost points from the source patch. For this to be the 
case, we need to introduce additional overlap points (in- 
dicated by colored boxes in Fig. |2| that lead to an over- 
lap of the evolved region. Effectively, this means that the 
equations are solved twice in the additional overlap re- 
gion, which introduces a small computational overhead. 

We note that quantities which are defined in the global 
Cartesian tensor basis such as the curvature evolution 
variables ([8| do not need to be transformed between lo- 
cal coordinates patches. In our present hydrodynamics 
scheme, however, the evolved conserved variables are de- 
fined in local coordinates. Hence, for inter-patch ghost 
zone interpolation, they must be transformed between 
local coordinate systems. Let us denote the local coor- 
dinates of source patch p as u^p^j, and the local coordi- 
nates of target patch q as coordinates u^^-j. The con- 
served density, which is a pseudo-scalar of tensor weight 
transforms as a scalar tensor density according to 



D. Cell-centered AMR and Refluxing 

We introduce cell-centered AMR in combination with 
a refluxing scheme at refinement level boundaries to en- 
sure conservation of rest mass and - in the absense of GR 
effects - also momentum and energy of the fluid |84[ I110| . 
Because gravity leads to sources and sinks for fluid mo- 
mentum and energy, these quantities are generally not 
conserved in curved spacetimes. This is reflected in the 
source terms of the fluid conservation laws ^ , which are 
zero only in flat space. The numerical fluxes in our flnite 
volume scheme between grid cells, however, must be con- 
served. Since we employ subcycling in time where flner 
grids take multiple small time steps for each coarse grid 
time step [81, the conservation properties of our flnite 
volume approach do not hold at mesh reflnement bound- 
aries without refluxing. 

In cell-centered AMR schemes, coarse cells are subdi- 
vided into multiple smaller cells, ensuring that coarse grid 
and flne grid cell faces align (see red line in the lower part 
of Fig. |4]). In contrast, the cell centers do not align. This 
is different from vertex-centered AMR schemes, where 
one aligns coarse and flne grid cell centers but not their 
faces (red line in the upper part of Fig. |4| . 

One may argue that vertex-centered schemes are more 
natural for wave-type equations such as the Einstein 
equations, which is why vertex-centered reflnement was 
originally implemented in the Carpet AMR driver. How- 
ever, refluxing requires cell-centered reflnement, and this 
comes with a certain added complexity that we describe 
below. 

a. Prolongation. Prolongation is the interpolation 
from coarse to flne-grid cells. In a vertex-centered scheme 
(and when assuming a reflnement factor of two), every 
second flne-grid point is aligned with a coarse-grid point. 
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7^ 



FIG. 4: Vertex-centered AMR (upper figure) versus cell- 
centered AMR (lower figure). In cell-centered AMR, two fine 
grid cell faces always coincide with a coarse grid cell face (red 
line). Thus, it becomes possible to sum up the two fine grid 
fluxes computed on cell faces to become one coarse grid cell 
flux. Cell-centered quantities, however, always need to be in- 
terpolated in the prolongation and restriction operation. In 
the vertex-centered case, every second grid point coincides 
with one coarse grid point. Thus, interpolation is not neces- 
sary for every point, and restriction becomes exact. 



and prolongation there corresponds to a copy. In be- 
tween coarse-grid points, one needs to interpolate. Cur- 
vature quantities are interpolated via a fifth-order La- 
grange polynomial. Hydrodynamics quantities are inter- 
polated via a second-order ENO interpolator [109] (also 
see Sec. IIC3) to avoid oscillations near discontinuities. 



In a cell-centered scheme, every fine-grid cell requires 
interpolation. We interpolate curvature quantities via a 
fourth-order Lagrange polynomial, and interpolate hy- 
drodynamics quantities via a second-order ENO interpo- 
lator. 

b. Restriction. Restriction transfers fine-grid infor- 
mation to the next coarser grid, after both have been 
evolved in time, and are aligned in time again. Differ- 
ent discretization errors will have led to slightly different 
results, and one overwrites the coarse-grid results by re- 
spective fine-grid results. For a vertex-centered scheme, 
this is straightforward, since each coaxse-grid point is 
aligned with a fine-grid point, and hence the variable on 
the fine-grid point can simply be copied. 

For cell-centered schemes, things are more complex, 
since restriction also requires interpolation. We interpo- 
late curvature quantities via a third-order Lagrange poly- 
nomial. Hydrodynamics quantities are averaged, corre- 
sponding to linear interpolation. This is a conservative 
operation, so that e.g. the mass in a coarse-grid cell is 
the sum of the masses in all contained fine-grid cells. 

The distinction between curvature and hydrodynamics 
quantities is crucial to achieving high accuracy. If one 
does not use higher-order operations for the curvature 
quantities, then the accuracy of the overall simulation 
is significantly reduced. On the other hand, one needs 
to employ a conservative interpolation scheme for the 
hydrodynamics quantities, but can accept a lower order 




At' < 
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FIG. 5: AMR time evolution, showing fluxes across cell faces, 
for both coarse (upper row) and fine cells (lower row). Time 
moves upwards; the fine grid (lower panel) takes multiple 
steps for each coarse grid step (upper panel). In Berger-Oliger 
AMR, the coarse and fine levels are evolved independently, 
and the sum of the fine grid fluxes crossing the green faces are 
not guaranteed to be equal to the coarse grid flux crossing the 
red face. At the end of a time step, the neighboring bold-faced 
coarse and flue cells may be in an inconsistent state, requiring 
refluxing to add a correction to the light blue coarse grid cell. 



of accuracy there. For restricting curvature quantities, 
we therefore use third-order polynomial interpolation. 

c. Refluxing. Refluxing is an algorithm to ensure 
conservation across mesh refinement boundaries [851 1110| . 
Since coarse and fine grids are evolved in time indepen- 
dently, it is not guaranteed that the fluxes leaving the flne 
grid are identical to those entering an abutting coarser 
grid (see Fig.js]). Refluxing integrates the coarse grid and 
flne grid fluxes across these faces, and then adjusts the 
coarse grid cell just outside the reflned region according 
to the flux difference. 

We outline the generic refluxing algorithm for a con- 
served quantity / in the steps below. 

1. We start with a flne grid level I + I and a coarse 
grid level I which are momentarily aligned in time, 
i.e. tl = 4^^! where i denotes the i-th step on the 
coarse level, and j denotes the j-th step on the flne 
grid. Due to subcycling in time, for any coarse- 
grid time step, there are twice as many flne-grid 
time steps, i.e. i — 2j. 

2. At the reflnement boundary (red line of Fig. |4] or 
red and green lines in Fig. [s]) , we store integrated 
coarse and flne grid flux registers P and l'"*"^ for 
some conserved quantity /. Due to the 2:1 mesh 
reflnement, there are four integrated flne grid flux 
registers for every integrated coarse grid flux regis- 
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ter. (Only two are visible in Fig. |4|. At t\ 
all registers are zero. 



3. Each refinement level is independently integrated 
forward in time until the two refinement levels are 
aligned in time again, i.e., until we have t\j^^ = 
^2j+2- During each integration step, the hydrody- 
namic evolution scheme computes fiuxes F for a 
quantity / located at all cell interfaces. At the re- 
finement boundary, we use the computed fine grid 
fluxes on the fine grid cell interfaces, and 

coarse grid fluxes on the coarse grid cell inter- 
faces to integrate coarse and flne grid flux registers 
forward in time, i.e., we independently integrate 



dtl 



i+i 



F 



i+i 



(23) 



at the refinement boundary. 



4. After restriction, when t\j^i ~ ^2^+21 '^^ 

and /' to compute a correction for conserved quan- 
tity /. The correction is obtained as follows. 

(a) The integrated fine grid flux register /'+^ is 
restricted to the coarse grid via 



^finc — '^^ 



(24) 



where TZ denotes the cell interface restriction 
operator. Note that since the flux registers 
are stored on cell faces, this operator is differ- 
ent from the operator used for the fluid state 
vector. 

(b) A correction Cj- for conserved quantity / on 
coarse grid level I is now obtained via 



C\ = ilL,-l')/A^x 



(25) 



where A' a; denotes the grid spacing of reflne- 
ment level I. 

5. The correction Cj- is added to the coarse grid cell 
on level I next to the reflnement boundary (blue 
cell in Fig. |5|, i.e. 



J corrected J f ' 



(26) 



This completes the refluxing operation. We repeat 
the steps 1-5 until the evolution is complete. 

The steps above are performed for any of the evolved 
conserved quantities D, 5'*, r, and Y^°'^. 

We note that the state thus obtained in the corrected 
coarse grid cells may be thermodynamically inconsistent, 
e.g., near the surface of a star, and may need to be pro- 
jected onto a self-consistent state. This is to be expected 
with our atmosphere treatment, as we discuss in Ap- 
pendix [O 



E. Time Integration and Multirate Runge-Kutta 

Schemes 



We carry out time integration using the Method of 
Lines (MoL) [5H] . MoL is based on a separate treatment 
of the spatial derivatives (the right-hand sides), and the 
time derivatives. This allows one to employ integration 
methods for ordinary differential equations (ODE) such 
as Runge-Kutta (RK) schemes for the time integration. 

We evolve the spacetime and hydrodynamic sector of 
our evolution system simultaneously using full matter- 
spacetime coupling. The coupling between the two sec- 
tors is achieved via source terms. The spacetime evolu- 
tion is sourced by the stress-energy tensor computed by 
the hydrodynamic sector. Vice versa, the hydrodynamic 
part contains additional source terms which are a result 
of the coupling to a curved spacetime metric. Written in 
simplified form, our system is given by 



dtg = F(g,q), 
dtH = G(g,q), 



(27) 
(28) 



where g denotes curvature evolution quantities, q de- 
notes fluid evolution quantities, and F and G denote the 
RHS functions. 

Traditionally, spacetime metric and hydrodynamic 
variables are evolved simultaneously using the same time 
integration scheme. A standard choice in our case is the 
classical fourth-order Runge-Kutta (RK4) method. The 
timestep is chosen such that the Courant-Friedrich-Lewy 
(CFL) factor, deflned as C = At/ Ax, becomes C = 0.4. 
The CFL factor is limited by the stability region of the 
numerical scheme, which in turn is limited by the speed 
of hght. 

We observe two important points in our simulations. 
First, the error in our numerical evolution is in most cases 



not dominated by the time integration (see Sec. Ill I. The 
choice of At is not guided by accuracy requirements, but 
rather by the restrictions imposed by the CFL condition. 
This is unfortunate since a larger timestep would speed 
up our simulation with only small negative impact on 
the accuracy. Second, we flnd that the CFL factor is 
largely determined by the spacetime evolution. In the 
Cowling approximation, i.e. when the spacetime sector is 
not evolved and held flxed at its initial setup, we typically 
can use more than twice as large CFL factors (up to 
C ~ 1) without encountering any numerical instabilities. 

Since our timestep is flxed, rather than enlarging the 
timestep At (and hence C), we switch to the classical 
second-order Runge-Kutta (RK2) method instead. This 
scheme has a smaller stability region by roughly a factor 
of two compared to RK4. Due to the less restrictive CFL 
factor for the fluid evolution compared to the curvature 
evolution, however, we can still use the same timestep 
as for the curvature evolution with RK4. The advantage 
of the RK2 schemes is that it require half as many RHS 
evaluations compared to RK4. The accuracy of RK2, 
however, is typically much lower than that of an RK4 
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TABLE II: Butcher tableau for an explicit multirate 
RK4/RK2 scheme. The right table (separated by the double 
vertical line) shows the coefficients bi (bottom line), d (first 
vertical column), and Oij for the classical RK4 scheme. The 
left table shows the corresponding RK2 coefficients evaluated 
at timesteps that coincide with RK4 timesteps. 
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1/2 




1/2 


1/2 




1/3 


1/6 


1/6 1/3 



integration. In practice, we find that the reduction in 
accuracy is not a severe limitation for most cases (see 
Sec.|mi. 

We therefore apply the RK2 integrator for the hydro- 
dynamic sector, while maintaining the RK4 integrator 
for the spacetime part. 

A scheme for coupling different parts of a system of 
equations with different RK integrators is given by mul- 
tirate RK schemes (e.g. [88| ^^). Here, we make the 
simple Ansatz of performing one RK2 intermediate RHS 
evaluation for two RK4 intermediate RHS evaluations. 
That is, the additional RK4 intermediate RHS evalua- 
tions simply use the results from the last intermediate 
RK2 step. 

To be more explicit, given the equation 



dty = f{t,v) , 



(29) 



where / corresponds to the RHS, we write a generic RK 
scheme according to 



(30) 



h = f{tn + CiM ,yn + l^t^a.jkj) . (31) 



The coefficients 6, , c, 



and ttij can be written in the stan- 



dard Butcher notation (see, e.g. |111| ). 

In our multirate scheme, we use two different sets of 
coefficients. The coefficients for the RK2 scheme are ar- 
ranged such that RHS evaluations coincide with RK4 
RHS evaluations. We list the corresponding multirate 
Butcher tableau in Table HTl 



F. Gravitational Wave Extraction 

GWs are extracted in the wave- extraction zone of our 
simulation. We define the wave-extraction zone as the 
region on the computational grid which is at sufficient 
distance from the gravitating source to avoid near-zone 
effects, and at the same time offers sufficient resolution 
to resolve the waves. Beyond the wave-extraction zone. 



we typically use radial stretching to gradually decrease 
the radial resolution up to a certain radius (e.g.. Fig. [T]). 

We use the techniques described in detail in |21j . 
Among those are (i) the standard slow-motion weak- 
field quadrupole formalism (see e.g., [151 [TBI [221 11121 - 
I114j ) which is purely based on the quadrupolar matter 
distribution and does not take into account any curva- 
ture effects, (ii) Regge-Wheeler-ZeriUi-Moncrief (RWZM) 
extraction based on gauge-invariant spherical perturba- 
tions about a fixed Schwarzschild background (see |115| 
for a review), (iii) Newman-Penrose extraction based on 
complex spin-weighted components of the Weyl tensor 
[731 11161 1117| , and (iv) Cauchy-characteristic extraction 
(CCE) making use of nonlinear nuUcone evolu- 

tions of the Einstein equations out to future null infin- 
ity (see [118, for a new high-order algorithm). The 
latter extraction technique is the only one capable of de- 
termining the gravitational radiation content unambigu- 
ously and without finite-radius and gauge errors [21-24]. 

The curvature-based techniques (ii)-(iv) require one or 
two integrations in time in order to compute the strain, 
which may lead to strong non-linear and unphysical arti- 
ficial drifts. This can be overcome by the fixed frequency 
integration (FFI) technique presented in |119| . FFI re- 
quires the choice of a cut-off frequency /o, which ideally 
must be below the physical frequency components con- 
tained in the signal. For instance, for a typical binary 
neutron star inspiral signal, /™ < mOoi.bitai/27r, where 
f^orbitai IS the initial orbital frequency, and ni is the as- 
sociated harmonic m-mode number. 

The energy and angular momentum that is lost due to 
the emission of GWs can be computed in terms of spin- 
weighted spherical harmonic coefficients of ^4 as derived 
in [12011121] . We use the expressions for the radiated en- 
ergy flux d-Eiad/c^i and angular momentum flux dJ^nd/dt 
in terms of the Weyl scalar ^4 from [121| . In the expres- 
sions for dEj-i^^/dt and dJ^i^^/dt, we evaluate the appear- 
ing time integrals of the harmonic modes using FFI with 
/q" = m/o for each given to- mode. In order to obtain the 
total radiated energy i?iad and angular momentum Jj-ad 
from their fluxes, respectively, we time integrate in the 
time domain^. 



1. Numerical Setup 

We report the numerical settings employed for the var- 
ious wave extraction techniques that are used in this 
work. Since we are not interested in the numerical 
convergence properties of the wave extraction methods 
themselves (this has been analyzed elsewhere, e.g. \21\ - 
[531 [Sni [73J [71] ) , we stick to flxed settings for all test cases 
and numerical resolutions considered in Sec. IIIII Guided 



^ FFI cannot be applied since the radiated fluxes are non- 
oscillatory. 
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by previous work f21U23| . we find that the numerical error 
in the wave extraction is negligible provided appropriate 
settings. 

The most involved GW extraction technique is CCE. In 
that method, we solve the Einstein equations along null 
hypersurfaces between a worldtube F and future null in- 
finity . The worldtube F is typically located at some 
radius i?r in the wave-extraction zone, and is simula- 
tion dependent [l^ (and references therein). Specific to 
the present work, the CCE grid consists of Nr = 301 
points along the radial direction. Each radial shell is 
discretized by two stereographic patches comprised of 
-^ang = 81 points per direction per patch. At the inner- 
boundary worldtube F, we use up to ^max — 8 harmonic 
modes for the decomposed Cauchy metric data. The 
metric data is decomposed on spheres with Nq — 120 
and = 240 points in 9 and direction, respec- 
tively. The compactification parameter^ r^t is set to the 
particular extraction radius for a given simulation, e.g. 
fvjt = IOOA/q. In all cases, the innermost radial com- 
pactified coordinate point is given by x^^ — 0.49. To- 
gether with an appropriate setting of r^t, this ensures 
that the worldtube F is located close to the first few ra- 
dial points on the characteristic grid. The timestep and 
extraction radius must be picked on a case by case ba- 
sis. The wave-extraction zone is alwaus located on the 
spherical "inflated-cube" grids. For the stellar collapse 



model A3B3G3 (Sec. IIIBl, the wave-extraction zone is 
located between radii 1000 M© < i?r < 25OOM0. For 
all remaining tests, the wave-extraction zone is located 
at 100 M0 < i?r < 25OM0. The wave-extraction output 
frequency is dictated by the timestep of the spherical 
"inflated-cube" grids. 

The remaining wave-extraction techniques are much 
simpler and only require single spheres at some finite ra- 
dius R. 

To project metric data from the 3D grid onto spheres, 
we use fourth-order Lagrange interpolation. 



TABLE III: Initial parameters and properties of the (per- 
turbed) TOV star used to construct the initial data. The 
density perturbation is only applied in the perturbed TOV 
test case. Units are in c = G = Mq = 1. 

Polytropic scale K 100 

Polytropic index F 2 

Central rest-mass density pc 1.28 x 10^^ 

ADM mass [Mq] Madm 1.4002 

Baryonic mass [A^q] Mb 1.5062 

Equatorial radius [Mq] ([km]) Re 9.586 (14.16) 



Density pert, mode 
Density pert, amplitude 
Monopole fundamental mode ]kHz 
First overtone ]kHz] 
Quadrupole fundamental mode ]kllz 
First overtone ]kHz] 



I 
A 
F 
Hi 

1 V 



2 

0.01 

1.458 

3.971 

1.586 

3.726 



which is part of the EinsteinToolkit. This framework 
defines mass and angular momentum in terms of partic- 
ular closed 2-surfaces, such as the apparent horizon. 

The spherical surface defining the apparent horizon 
shape uses Ng = 41 points along the 0-direction and 
Na, = 80 points along the 0-direction. 



III. RESULTS 

We revisit a number of "benchmark" problems com- 
monly found in the literature: an isolated perturbed and 
unperturbed neutron star, a rotating core collapse model, 
a collapsing neutron star to a black hole, and a binary 
neutron star coalescence. Basic code tests such as shock 
tubes can be found in the Appendix. We describe our 
analysis in more detail in corresponding sections below. 



Isolated Neutron Star 



G. Horizon Finding and Hydrodynamic Excision at 
the Puncture 

To track the appearance and shape of an apparent hori- 
zon, we use AHFinderDirect [IS] which is part of the 
EinsteinToolkit [83J. As soon as an apparent horizon 
is found during an evolution, we excise the fluid variables 
within a fraction of the radius of the apparent horizon 
and set them to their corresponding atmosphere values. 
We get stable evolutions when excising about 85% of the 
interior of the apparent horizon volume. 

In order to compute angular momentum Jah and mass 
A/ah of a black hole, we use the isolated / dynamical hori- 
zon framework provided by QuasiLocalMeasures |122| . 



^ See 1231 for a description of CCE relevant parameters. 



We investigate convergence and accuracy of an iso- 
lated unperturbed neutron star and an isolated perturbed 
neutron star using full GR matter-spacetime coupling in 
three spatial dimensions. The neutron stars are given by 
the solution of the Tolman-Oppenheimer-Volkoff (TOV) 
equations fnHlIT^ . 

This test aims at showing the correctness of our cell- 
centered AMR scheme, and enhanced PPM reconstruc- 
tion. 



1. Initial Conditions and Equation of State 

We use a polytropic equation of state P — Kp^ with 
scale K = 100 and index F = 2 in the initial data con- 
struction. Although this choice does not represent a real- 
istic choice for real neutron stars, these parameters have 
been used in previous work (e.g. |361 1125| ). and can be 
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FIG. 6: Unperturbed TOV star: normalized central density 
Pc{t)/pc{t = 0) — 1 on the three resolutions rO, rl, and r2 
(top panel), difference in normalized central density between 
low and medium resolutions, and medium and high resolu- 
tion (center panel), and the L2-norm of the Hamiltonian con- 
straint ||-ff||2 on all three resolutions (bottom panel). As the 
resolution is increased, the amplitude of the central density 
oscillations, the offset, and the slope decrease as expected. 
The differences in resolutions of the central density are scaled 
for second-order convergence. The L2-norms of the Hamilto- 
nian constraint ||-ff||2 are scaled for first-order convergence. 
The resolution study is performed using cell-centered AMR 
and ePPM. 




o 



-12 



ePPM, cc 

■ oPPM, vc 



oPPM, cc 
ePPM, cc, multirate 




ePPM, cc 

oPPM, vc 



oPPM, cc 
ePPM, cc, multirate 



0.0 



1.5 3.0 



4.5 6.0 
t [ms] 



7.5 9.0 



FIG. 7: Unperturbed TOV star: the L2-norm of the Hamilto- 
nian constraint ||-ff||2 (upper panel), and conservation of total 
baryonic mass Mb (lower panel) for different numerical se- 
tups. We compare vertex-centered (vc) with cell-centered (cc) 
AMR using oPPM and/or ePPM. In addition, we also show a 
simulation with "ePPM, cc" using multirate time integration. 
\\H\\2 is strongly effected by the choice of numerical scheme, 
while Mb is essentially uneffected. The setup "ePPM, cc" 
performs best, while "oPPM, cc" performs worst. The stan- 
dard scheme "vc, oPPM" used in other codes (e.g. [83! !S2]) 
is slightly worse than the new scheme "ePPM, cc". Multirate 
time integration leads to nearly identical results. 
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FIG. 8: Perturbed TOV star: the top panel shows the "+" po- 
larization of the GW strain Z)/i+,e as emitted in the equatorial 
plane and rescaled by distance D for the three resolutions rO, 
rl, and r2. The waveforms are computed with CCE. In the 
panel below, we show the differences in GW strain between rO 
and rl, and rl and r2, where the latter is rescaled for second- 
order convergence. In the third panel from the top, we show 
the absolute central density evolution pc{t) for the three res- 
olutions. Below, we show the differences in central density 
scaled for second-order convergence. In the bottom panel, 
we show the L2-norms of the Hamiltonian constraint |j-?/|l2. 
Since the initial data for the perturbed case are not constraint 
satisfying, the constraints do not exhibit clean convergence. 
The convergence study is performed using cell-centered AMR 
and ePPM. 



used as code verification. During evolution, we use an 
ideal fluid F-law equation of state with F = 2. The 
key parameters are given in Table |III| The initial data 
are generated via Hachisu's self-consistent field method 
[1261 1127j which requires as input the central density pc 
of the star, and a polar-to-equatorial axes ratio between 
and 1 to define rotation. In the present case, we set 
Pc — 1.28 X IQ-'^ and use an axes ratio of 1 (no rota- 
tion). In the case of the perturbed TOV star, we perturb 
the star by a spherical harmonic {£,m) = (2,0) density 
perturbation of amplitude A = 0.01. 



2. Numerical Setup 

The grid is similar to the one depicted in Fig. [T] ex- 
cept that here, we have just one refinement region. The 
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FIG. 9: Perturbed TOV star: the impact of different numer- 
ical settings on the I/2-norm of the Hamiltonian constraints 
\\H\\2 (upper panel), and on the conservation of baryonic mass 
Mb (lower panel) . The setup using vertex-centered (vc) AMR 
and oPPM (blue dashed curve) leads to larger constraint vi- 
olations than the setup using cell-centered (cc) AMR and 
ePPM. Multirate time integration does not change the ac- 
curacy of the results. In all cases, Mb is nearly equally well 
conserved. 
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FIG. 10: Perturbed TOV star: Power spectrum of poo, P20 
and /i+,e (individually scaled for better visibility), and the 
first few fundamental neutron star oscillation modes (vertical 
lines) computed in |36] . 



fine grid spacing is Ax = 0.2 Mq for the low resolu- 
tion (rO), Ax = 0.125 for the medium resolution 
(rl), and Aa; — 0.1 Mq for the high resolution simu- 
lation (7-2). The fine grid extends to i? = 11 Mq and 
encompasses the entire star. The inter-patch bound- 
ary between central Cartesian patch and outer spher- 
ical grid is located at i?s = 65 M©. We use 15, 24 
and 30 cells per angular direction per spherical patch 
for the low, medium and high resolutions, respectively. 
The radial resolution is chosen based on the Cartesian 
coarse grid resolution Ar = 1.6 Mq, Ar = 1.0 Mq, 
and Ar — 0.8 M©, for low, medium, and high resolu- 
tions, respectively. We use radial stretching outside the 
wave extraction zone to efficiently extend the computa- 



tional domain so that the outer boundary is causally dis- 
connected from wave-extraction zone and interior evo- 
lution. Accordingly, we stretch the radial resolution to 
Ar = 6.4 Af©, Ar = 4.0 M©, and Ar = 3.2 M© for low, 
medium, and high resolution simulations, respectively, in 
the region between radii i?i = lOOM© and i?2 — 800M©. 
The outer boundary is located at i?B — 3500 M© . 



3. Discussion 

a. Unperturbed TOV star. We first consider a sin- 
gle isolated non-rotating, unperturbed TOV star with 
parameters reported in Table In the top panel of 
Fig. [6] we show the normalized central density evolution 
Pc{t)/pc{t = 0) as a function of time on the three reso- 
lutions rO, rl, and r2, using our new cell-centered AMR 
and enhanced PPM scheme. In an ideal setting, the cen- 
tral density evolution should be constant as a function 
of time since the TOV solution represents a static fiuid 
configuration. Numerical errors induced by interpolation 
from the initial data solver grid onto the evolution grid, 
however, lead to an artificial excitation of the star, and, 
hence, to non-trivial central density oscillations, which 
must converge to zero as the resolution is increased. Due 
to the interpolation of the fiuid initial data onto the evo- 
lution grid, we observe a large initial spike and an over- 
all offset in the density oscillations. We additionally see 
an overall non-zero slope in the central density evolution 
caused by numerical errors during evolution. As the res- 
olution is increased, we consistently observe that the am- 
plitudes of the oscillations decrease, the offset becomes 
smaller, and the overall slope is reduced. In the center 
panel, we show the difference in normalized central den- 
sity pc{t)/pc{t = 0) between resolutions rO and rl, and 
rl and r2. We perform a three-level convergence test by 
computing the ratio of the differences in a given quantity 
F between the three resolutions, 
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(32) 



The ratio C defines the measured convergence rate of the 
solution (e.g. |103j ). Given three resolutions with spac- 
ing Axiow, Axnicdium, and Aa;high, the theoretical conver- 
gence rate for a particular order of convergence p can be 
computed via 



|A< 



medium 



p 

'^high 



(33) 



medium i 



Given our numerical resolutions, according to (33 1, we 
expect that the difference between medium and high res- 
olution, rl and r2, decreases by a factor of C = 4.33 for 
second-order convergence compared with the difference 
between medium and low resolution, rl and rO. 

In the bottom panel of Fig. [6] we show the time evo- 
lutions of the L2-norm of the Hamiltonian constraint 
||i?(i)||2 (16) for the three resolutions rO, rl, and r2. 
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As the resolution is increased, the error drops consistent 
with first-order convergence, since the rescaled medium 
and high resolution curves are on top of each other. We 
note that while the fluid body itself is smooth, the sur- 
face of the star is non-smooth, hence inducing a dominant 
flrst-order error (compare Fig. 25 1. 

In the top panel of Fig. [7j we show the L2-norm of the 
Hamiltonian constraint II2 of a static TOV star using 
vertex-centered (vc) AMR, and cell-centered (cc) AMR. 
Both AMR setups are run with the oPPM and ePPM 
reconstruction method. In addition, we also perform a 
simulation using cell-centered AMR and ePPM recon- 
struction with multirate time integration. We observe 
that the setup "ePPM, cc" exhibits the lowest constraint 
violations. The setup "ePPM, cc, multirate" is right on 
top of the red curve, hence indicating comparable accu- 
racy. The setup "vc, oPPM", which is the setup used in 
previous work (e.g. [U El EU |83l [92] ) yields slightly less 
accurate evolution. Finally, the setup "cc, oPPM" yields 
significantly reduced accuracy compared to all other se- 
tups. This is mainly due to the oPPM scheme, which is 
known to reduce the order of accuracy at smooth max- 



ima to first order (see Appendix \S\ Fig. 25 ) . This effect 



is not seen in the vertex-centered setup "vc, oPPM", since 
the central density is exactly located on a grid point. 

In the bottom panel of Fig. |7] we show conservation of 
mass for the considered numerical setups. In all cases, the 
total mass loss is on the order of 10^^ over the course of 
the evolution. Since the AMR boundaries are all located 
in the vacuum region outside the star, refluxing at AMR 
boundaries is not relevant. The mass loss is entirely due 
to interaction with the artificial low-density atmosphere 
in the vacuum region (see also Appendix [C| . 

b. Perturbed TOV star. As a second test, we apply 
an initial {£ = 2, m) = (2, 0) density perturbation with 
amplitude A = 0.01 onto the same TOV star considered 
above. A more complete study of this configuration in- 
cluding variations on perturbation parameters has been 
performed in |361 1125] . Numerical grids and setups are 
identical to those of the static TOV star, and we perform 
the same analysis as above. In addition, we also analyze 
the non-trivial {£,m) = (2,0) mode of the GW signal 
that is induced by fundamental mode oscillations. In the 
upper panel of Fig. [8] we plot the "+" polarization of 
the GW signal Dh+ f. as emitted in the equatorial plane 
from the three resolutions rO, rl, and r2. Since only the 
{£, m) = (2, 0) mode is excited, the entire wave signal can 
be written as 



Dh+,, ^ Dhf _2Y2o{0 = 



= 0). 



(34) 



Here, D is the distance from the source. We compute 
with CCE and u se an FFI cut-off frequency of /q = 812 
Hz (see Sec. II F I . We also show the differences of the GW 
strain between low and medium, and medium and high 
resolutions, where the latter is scaled for second-order 
convergence. In addition, we show the central density 
evolution Pc(t) for the three resolutions which converge. 
Similar to the above, we plot the differences between low 



and medium, and medium and high resolutions scaled for 
second-order convergence. We also show the L2-norm of 
the Hamiltonian constraints \\H\\2 of the three resolu- 
tions. Since the initial data solver does not take into 
account the effects of the perturbation onto the initial 
spacetime metric, the constraints do not converge ini- 
tially, and only slowly converge at later times. In the 
present plot, we have not used any rescaling. We note, 
however, that the slopes of the medium and high reso- 
lutions are slightly smaller than for the low resolution 
case. 

When comparing the strain Dh'^^ as computed with 

CCE to the strain Dh'^^ as computed from the RWZM 
formalism, we generally find that the strain computed via 
the RWZM formalism is prone to numerical noise. In ad- 
dition, we find that the finite-radius error and gauge error 
inherent in the waveform obtained from RWZM master 
functions at radii R = 100 Mq and R = 250 Mq is on the 
order of 10%. A similar behavior applies to the strain 
Dh}^^ as extracted via the NP formalism at a finite ra- 
dius. 

Finally, we also check that the correct fundamental 
oscillation modes are excited. In Fig. [lO] we compare 
the frequency spectrum of the density p and the strain 
-D/i+^e to the eigenmodes found in [36 . In order to com- 
pute the spectrum of p, we first project p from the 3D 
grid onto spherical shells inside the star, and then de- 
compose in terms of spherical harmonics. The vertical 
lines in Fig. [TO] correspond to the fundamental monopole 
mode F and its first overtone Hi, and the fundamen- 
tal quadrupole mode ^/ and its first overtone ^pi. As 
expected, the spectrum of the strain Dh^ f, and the 
{£, m) = (2, 0) mode of the density p20 both peak at the 
correct quadrupole eigenmode frequencies. Likewise, the 
spectrum of the {i, m) = (0,0) density mode correctly 
peaks at the monopole eigenmode frequencies. 



B. Rotating Stellar Collapse 

We investigate convergence and accuracy of the bench- 
mark rotating stellar collapse model A3B3G3, which has 
been previously considered in the literature |28| [25] . This 
tests the ability of the code to simulate the collapse of 
a rapidly differentially spinning iron core in full 3D with 
causally disconnected outer boundaries, albeit with sim- 
plified microphysics. We show that due to larger wave 
extraction radii, the waveforms extracted via curvature- 
based methods such as CCE are more accurate than what 
has been computed before [21^. 



1. Initial Data and Equation of State 

For the purpose of this test, we employ a hybrid 
equation of state |29| [30] I128| that combines a 2-piece 
piecewise polytropic pressure Pp with a thermal com- 
ponent Pth, i.e., P = Pp + Pth- To model the stiffen- 
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FIG. 11: Stellar collapse: The GW strain Dh+^e extracted 
via the quadrupole formula (upper panel), the central density 
Pc (third panel from the top), and the I/2-norm of Hamilto- 
nian constraint ||-ff||2, all on the three resolutions rO, rl, and 
r2. The panels directly below the top panel and the third 
panel from the top show the difference in strain and central 
density between low and medium resolutions, and medium 
and high resolutions. The differences are scaled for second- 
order convergence. The L2-norm of Hamiltonian constraint 
is scaled for first-order convergence. Before core bounce, the 
constraint exhibits second order convergence. After shock for- 
mation, the convergence rate is reduced to first order. The 
convergence study is performed using cell-centered AMR and 
ePPM. 



ing of the equation of state at nuclear density 
2 X lO^^'gcm"^, we assume that the polytropic index T 
jumps from Fi below nuclear density to F2 above. The 
equation of state parameters are given in Table |IV| 

The initial data are constructed from n = 3 (Fijni — 
Fi — 4/3) polytropes in rotational equilibrium gener- 
ated via Hachisu's self-consistent field method |126| 1127] 
which not only provides fluid, but also spacetime cur- 
vature initial data. While being set up as marginally 
stable polytropes with Fi.ini = 4/3, during evolution, 
the initial sub-nuclear polytropic index Fi is reduced to 
Fi < Fi i„i to accelerate collapse. Following previous 
studies |16l we use r2 — 2.5 in the super-nuclear 
regime. 

In the present test, we revisit model A3B3G3 from 
[28, '29j. This conflguration uses Fi = 1.31. It is strongly 
differentially rotating, with its initial central angular ve- 
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FIG. 12: Comparison of vertex-centered (vc) AMR with 
oPPM versus cell-centered (cc) AMR with ePPM for stellar 
collapse model A3B3G3. We show the central density pc (up- 
per panel), tthe L2-norm of the Hamiltonian constraint ||-ff||2 
(middle panel), and conservation of total baryonic mass Mb 
(bottom panel). Due to refluxing in the cell-centered case, 
the mass is almost perfectly conserved, while in the vertex- 
centered case, the mass is rapidly growing (bottom panel). 
Due to ePPM, the constraints in the cell-centered case ex- 
hibit almost no growth after core bounce, while in the vertex- 
centered case with oPPM the constraints are clearly growing 
(lower panel). The results are not changed when multirate 
time integration is used. The comparison is done using base- 
line resolution rl. 



TABLE IV: Initial parameters and properties of the rotating 
stellar collapse model A3B3G3. Units are in c = G = Mq = 
1. 



Polytropic scale 


K 


0.4640517 


Initial polytropic index 


ri,ini 


1.3 


Evolved polytropic index 1 


Ti 


1.31 


Evolved polytropic index 2 


r2 


2.5 


Thermal polytropic index 


Pth 


1.5 


Central rest-mass density 


pc 


1.6193 X 10^ 


Axes ratio 




0.93 


Degree of differential rotation [km] 


A 


500 


Rotational / binding energy [%[ 


T/\W\ 


0.9 


Equatorial radius [Mq] 


Re 


1.0661 X 10^ 


Baryonic mass [Mq] 


Mb 


1.4596 


ADM mass [Mq] 


Madm 


1.4596 


ADM ang. mom. [Mq] 


Jadm 


2.4316 


Spin 


a 


1.1413 



locity dropping by a factor of two over A = 500 km. This, 
in combination with T/|Ty| = 0.9%, leads to rapid rota- 
tion in the inner core, resulting in a very strong GW 
signal at core bounce and dynamics that are significantly 
affected by centrifugal effects. It produces a "Type-I" 
GW signal with a centrifugally-widened broad peak at 
core bounce [28i |29j . 
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FIG. 13: The GW strain extracted Dh+,e from the rotat- 
ing stellar collapse model A3B3G3 (upper panel). We show 
the strain extracted via CCE from different worldtube loca- 
tions Rr = 1000 M©, _Rr = 1500 Mq, and Rr = 2500 M©, 
as well as the strain computed via the quadrupole formula. 
Larger CCE worldtube radii permit lower FFI cut-off frequen- 
cies without introducing unphysical drifts in the GW strain. 
All waveforms extracted via CCE are in good agreement to 
within a few percent with the waveform computed via the 
quadrupole formula. The lower panel shows the differences in 
strain amplitude of the inner extraction radii to the outermost 
extraction radius. The differences converge as the extraction 
radius is increased. The comparison is done using baseline 
resolution rl. 



2. Numerical Setup 

We use five refinement levels located at the 
center of the domain. The refinement boxes 
of each level have a half-width of R^y = 
[192M0 , 144M0 , 98M0 , 4OM0 , 12M0] , respectively. 
The coarsest level is comprised of cubed-sphere 
multipatch grids (Fig. [TJ. The inner radius of the 
spherical grids is i?g = 384Af0, and the outer bound- 
ary is i?B — I6OOOM0. Initially, only the coarsest 
level is active. Additional levels are progressively 
added as the central density increases during col- 
lapse. The initial stellar radius of model A3B3G3 is 
Re — 1066. IM0 — 1574.84 km in the equatorial plane. 
Thus, the inter-patch boundaries thread the star in 
this particular setup. The finest refinement level is 
picked such that the protoneutron star is fully contained 
on that level. The GW extraction zone extends to a 
radius oi R = 25OOM0. Beyond that radius, we apply 
radial stretching up to a radius R — 6000 Af0. In this 
stretching region, the radial grid spacing is increased by 
a factor of 16, and the resolution becomes too coarse for 
reliable wave extraction. 

For our baseline resolution (denoted by rl), we pick a 
radial grid spacing of Ar = 8.OM0 on the non-stretched 



spherical inflated-cube grids, and a Cartesian resolution 
of Ax = 8.OM0 on the central Cartesian patch. Given 
our five refinement levels above, this results in a reso- 
lution of O.25M0 — 369.3 m for the protoneutron star. 
The angular resolution of the cubed-sphere grids is set 
to A^ang — 30 cells per patch and direction. This makes 
a total of A'ang, total — 120 points across the equatorial 
plane. 

In addition to our baseline resolution rl, we also con- 
sider a low resolution run rO, and a high resolution 
run r2 to check for convergence. Resolution rO uses 
Ar = Ax = 9.6Af0 and iVa„g = 24 (20% lower), and 
resolution r2 uses Ar = Ax = 6.4M0 and iVang = 36 
(20% higher). 

In all considered cases, we set the damping coefficient 
of the T-driver gauge condition to rj = 1/2. Dissipation 
is set to ediss = 0.1 on the fine levels, and e^iss — 0.01 
on the multipatch grid. The atmosphere level is set to 
be 10^^° of the central density, and we damp the stress- 
energy tensor in the atmosphere starting using ([7| with 
Ra = 13OOA/0 and i?i = 14OOM0. 



3. Discussion 

In Fig. |11| we show convergence of the plus polariza- 
tion of the GW strain Dh^ ^ measured in the equatorial 
plane, the central density pc, and the X2-norm of the 
Hamiltonian constraint ||-ff||2- The GW strain is com- 
puted using the quadrupole formula, though a similar 
analysis and result applies to all extraction methods. All 
three quantities are shown for the three resolutions rO, 
rl, and r2, using multipatches, cell-c entered AMR, re- 
fiuxing, and enhanced PPM (see Sec. Ill B 2 1. We align 



the results from all three resolutions at the time when the 
central density pc reaches its maximum at core bounce. 
We observe first order convergence in \\H\\2 after core 
bounce. In the prebounce phase, \\H\\2 exhibits second- 
order convergence. This behavior is expected since the 
numerical scheme reduces to first order at the shock front 
after bounce where the error are greatest. 

In Fig. [TT] we also show the absolute difference of the 
GW strain D/i+.e and the central density pc between low 
(rO) and medium (rl) resolutions, and medium and high 
(r2) resolutions. The convergence behavior of the two 
quantities is less clean than what can be observed for 
the Hamiltonian constraint due to their oscillatory na- 
ture. The convergence is between the expected first and 
second-order accuracy. 

In Fig. [12] we compare vertex-centered AMR with 
original PPM reconstruction versus cell-centered AMR 
with refiuxing and enhanced PPM. In addition, we show 
the behavior of the latter case when multirate RK time 
evolution is applied. As is clear from the bottom 
two panels, the cell-centered scheme with refiuxing and 
enhanced PPM ("cc, ePPM") outperforms the vertex- 
centered scheme with original PPM ("vc, oPPM"). While 
in the cell-centered case, ||i?||2 essentially remains con- 
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stant after core bounce, it clearly grows in the vertex- 
centered case. Even worse, the vertex-centered case ex- 
hibits a rapid growth in total baryonic mass after core 
bounce. The evolution with multirate RK performs 
equally well as the "cc, ePPM" setup, which uses standard 
RK4 time integration. The multirate setup offers a speed 
up of ~ 20% for the current test problem. The speed- 
up can be significantly larger when full microphysics and 
neutrino transport is employed (e.g. [OOj). 

In Fig. |13| we revisit our study of extracting gravi- 
tational radiation using curvature-based methods |21| . 
In |21| . we found a radial dependence of the accuracy 
of the curvature-based extraction methods. This study 
made use of purely Cartesian simulation domains, and 
was thus limited in terms of possible domain sizes and 
extraction radii. The maximum extraction radius was 
limited to i? = 1000 Af©. This is still fairly close and 
means that the waveforms are extracted well inside the 
star. Our curvature-based extraction methods, however, 
assume vacuum, i.e. a vanishing stress-energy tensor at 
the extraction location. In [21 , we thus conjectured 
that increased extraction radii that are located outside 
the star would further improve the accuracy of the ex- 
tracted waveforms. Given our new multipatch setup, we 
can confirm this conjecture. We have placed three ex- 
traction radii at i? = [lOOOM©, ISOOM©, 25OOM0] in a 
region with constant radial spacing Ar = 8.OM0 where 
the radial direction is not yet stretched. The upper panel 
of Fig. [13] shows the "+" polarization of the GW strain 
Dh^^e measured in the equatorial plain extracted via 
CCE. As a comparison, in the same panel, we also show 
Z?ft.+.e computed via the quadrupole formula. We ap- 
ply FFI to compute the strain Dh from ^4 extracted 
with CCE (see Sec. IIFl. In [3T], we conjectured that 
the low cut-off frequency that must be picked for FFI can 
be reduced as the extraction radius is increased. Here, 
we confirm that this is indeed the case. While extrac- 
tion radius R = IOOOMq requires a low cut-off frequency 
/o = 100 Hz which is well inside the LIGO sensitivity 
band, we find that at radius R = ISOOM© we can get 
away with fo — 60 Hz. At radius R = 25OOM0, we can 
further reduce this to /o = 30 Hz without introducing ar- 
tificial non-linear drifts in the strain. In the bottom panel 
of Fig. [13] we show the difference in GW amplitude of the 
waveforms computed from the inner extraction radii to 
the waveform computed from the outer most extraction 
radius. We confirm that as the extraction radius is in- 
creased, the differences further decrease similar to what 
has been found in PI] • 

The waveform computed via the quadrupole formula 
does not suffer from amplification of low frequency errors 
[21]. We observe that the waveforms extracted via CCE 
at larger radius and decreased /q more closely resemble 
the monotonically rising signal in the prebounce phase 
that the waveform computed via the quadrupole formula 
exhibits. Overall, in accordance with [21], we still mea- 
sure the same deviations between GW amplitudes com- 
puted from CCE and the quadrupole formula to within a 



few percent at core bounce. This is not surprising, since 
the error in CCE due to different worldtube extraction lo- 
cations is much smaller than the observed deviation from 
the waveform extracted via the quadrupole formula. 

Finally, we note that we have also computed the GW 
strain via the RWZM formalism (not shown). In our 
previous more detailed study on GW extraction in the 
context of rotating stellar collapse [H], we found that 
the RWZM formalism leads to waveforms which are con- 
taminated by high frequency noise. Unfortunately, in 
the current study, which allows us to use larger extrac- 
tion radii than R = IOOOM0, we find that the systematic 
high-frequency noise inherent in the RWZM waveforms 
is not reduced, but instead even increases with increased 
extraction radius. As already conjectured in [21], this is 
most likely due to the perturbative manner the waves are 
extracted from the spacetime in the RWZM formalism. 
In this formalism, the spherical background geometry is 
projected out, which can result in very small values for 
the aspherical perturbation coefficients that are prone to 
numerical noise and cancellation effects. At larger radii, 
the aspherical perturbations are even smaller since they 
fall of as 1/r, and thus are harder to capture accurately. 
The RWZM approach may therefore be less suited for the 
extraction of the generally weak GW signals emitted in 
core collapse. 



C. Neutron Star Collapse 

Three-dimensional collapse of an isolated neutron star 
to a black hole is a valuable test of accuracy and con- 
vergence of our code for black hole formation in massive 
stars. We consider the uniformly rapidly rotating model 
_D4 previously studied in |51|92J as a benchmark problem. 
Apart from showing convergence and consistency with 
previous results, we improve the simulations by causally 
disconnecting the outer boundary from the interior evo- 
lution and the wave-extraction zone. We show that cell- 
centered AMR with refluxing leads to better conservation 
of mass than vertex-centered AMR. We also employ CCE 
for GW extraction. 



1. Initial Data and Equation of State 

The initial condition is given by a stable relativistic 
polytrope. Specifically, we use a polytrope P — Kp^ 
with r = 2 and Ki^i = 100 in the initial data construc- 
tion. The initial data are generated via Hachisu's self- 
consistent field method |126[ 1127] . The central density 
is set to pc = 3.116 x 10"^ = 1.924 x lO^^gcm-^. We 
use an axes ratio of 0.65, which results in /3 = T/|VK| = 
7.6796 X 10^^ corresponding to a dimensionless spin of 
a = J/Af^ — 0.54354. In order to induce the gravita- 
tional collapse, we introduce an artificial pressure deple- 
tion of 2% by setting 7^ = 98 at the onset of the evolution. 
During evolution, we use an ideal fluid F-law equation of 
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FIG. 14: Rotating neutron star collapse: convergence analysis 
of the "+" polarization of the GW strain -D/i+,e as emitted in 
the equatorial plane and extracted via CCE (top two panels), 
central density pc evolution (next two panels), and L2-norm 
of the Hamiltonian constraint \\H\\2 (bottom panel). The dif- 
ferences in Dh+^e and pc between medium and high resolution 
are scaled for second-order convergence. At t — tBH = 0, the 
density drops to zero due to hydrodynamic excision within the 
horizon. The L2-norm of the Hamiltonian constraint (bot- 
tom panel) does not converge initially due to numerical ar- 
tifacts from the initial data solver, however, later converges 
at second-order during black hole formation t — ten ~ and 
black hole ring-down t — tsn > 0. The convergence study is 
performed using cell-centered AMR with ePPM. 



state with F — 2. The initial parameters and properties 
of the test case are summarized in Table IVl 



2. Numerical Setup 

The GW extraction is carried out on the cubed-sphere 
grid setup shown in Fig. [T] We pick the radius of the 
outer boundary such that the wave-extraction zone and 
the interior evolution are causally disconnected from the 
outer boundary, which we set to i?B — SOOMq. 

For our baseline grid setup rl, we make use of a radial 
and Cartesian resolution of Ar — Ax — 1.28M0 and 
N., 
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25 cells per patch and per angular direction. The 
boundary between central Cartesian and cubed-sphere 
grids is located at i?s = 65Af0. The radial coordinate 
spacing is increased from Ar to 2Ar in the region between 
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FIG. 15: Rotating neutron star collapse: we compare vertex- 
centered (vc) AMR and oPPM reconstruction with cell- 
centered (cc) AMR and ePPM reconstruction. The latter 
setup is also shown using multirate RK time integration. The 
top panel compares the central density evolution profile Pc{t). 
The center panel compares the evolution of the L2-norm of 
the Hamiltonian constraint ||-ff||2. The bottom panel com- 
pares the conservation of baryonic mass AIb- The setup "vc, 
oPPM" produces slightly larger violations in the Hamiltonian 
constraints, especially in the late collapse phase shortly before 
the black hole forms. Due to refluxing, the cell-centered case 
exhibits much better conservation of baryonic mass. Multi- 
rate RK time integration does not lead to different results. 
The comparison is done using baseline resolution rl. 



R = 25GM0 and R = 6OOM0. 

We employ five additional levels of AMR with half- 
widths R,i = [3OM0,18M0,11M0,5M0,3M0] located 
at the center of the Cartesian domain. With an ini- 
tial radius of i?NS ~ IOM0 along the equatorial plane, 
this means that the finest two levels thread through the 
neutron star. These two levels are required to resolve 
the black hole formed in the collapse. For our base- 
line resolution rl, we therefore have a grid spacing of 
Aa; = O.I6M0 = 0.24 km on the third finest level en- 
compassing the entire neutron star, and a resolution of 
Aa; = O.O4M0 — 0.06 km on the finest level containing 
the black hole. 

In addition to rl, we also use a low resolution rO with a 
coarse grid spacing of Ar = Aa; — 1.6Af0 and A^ang = 20 
cells per patch and per angular direction, and a high 
resolution setup r2 with a coarse grid spacing of Ar = 
Aa; = l.O24Af0 and A^ang = 31 cells per patch and per 
angular direction. 

We set the damping coefficient of the F-driver gauge 
condition to 77 = 1/2, and exponentially damp 77 to zero 
starting from radius i?^ = 65Af0 . 

The artificial low-density atmosphere is 10~* of initial 
central density. We also perform a simulation with an 
atmosphere density lO^^*' of the central density, however. 
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FIG. 16: Rotating neutron star collapse: we show the total 
ADM mass Madm (top panel; red dashed line) and the mass 
of the apparent horizon Mah plus energy radiated in GWs 
-Brad as a function of time (blue straight line). The total 
ADM angular momentum Jadm of the spacetime (red dashed 
line), and the angular momentum J ah as measured on the 
apparent horizon (blue straight line) is shown in the bottom 
panel. The inset plots show a close-up of the time evolution 
of Mah + Srad and J ah- As all matter becomes trapped 
in the event horizon, both, Mah + Emd and J ah, quickly 
asymptote to the conserved ADM values of the spacetime. 
Due to systematic (atmosphere) and numerical errors, the 
asymptoted values do not agree with the initial ADM values. 
Note that the mass radiated in GWs is negligible compared to 
the total mass of the black hole and thus barely contributes to 
Mah + -Erad • No angular momentum is radiated in GWs. The 
results are shown for resolution r2 using cell-centered AMR 
with ePPM. 



we find only negligible differences in the accuracy of our 
results. 



3. Discussion 

Following initial pressure depletion, the uniformly ro- 
tating polytrope collapses. During collapse, the central 
density pc increases until time t — t-BH — 0, the time when 
an apparent horizon, and thus a black hole forms. After 
formation of the horizon, the matter inside the horizon is 
excised from the grid, and the remaining exterior matter 
is rapidly dragged into the nascent black hole, leaving 
behind the artificial low-density atmosphere. Upon for- 
mation, the black hole is highly excited and radiates GWs 
until it settles to a Kerr state. This produces a character- 
istic ring-down GW signal with a particular quasi-normal 
mode frequency which depends only on mass and spin of 
the black hole. 

we show the emitted GW signal D/i+.e 



In Fig. 14 



and the evolution of the central density Pc for the three 
resolutions rO, rl and r2, respectively. The simulations 
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FIG. 17: Rotating neutron star collapse: power spectral den- 
sity of "+" polarization of GW strain D/i+.c as emitted in the 
equatorial plane and extracted via CCE. The blue straight 
line is the spectrum of the entire waveform, while the green 
dashed line is the spectrum of the ring-down signal. The red 
vertical line denotes the {£, m) = (2, 0) prograde fundamental 
{N — 0) quasi-normal mode frequency /qnm ~ 6.68 kHz of 
a spinning black hole of mass M = 1.8602 Mq and dimen- 
sionless spin a = 0.5435 as computed in |129| . Mass and 
spin of the nascent black hole are determined on its apparent 
horizon using the isolated horizon framework. The analysis 
is done using baseline resolution rl with cell-centered AMR 
with ePPM. 



are performed using cell-centered AMR, refiuxing, and 
ePPM reconstruction. The GW signal is extracted us- 
ing CCE and we use FFI with a cut-off frequency of 
/o = 1 kHz to obtain Dh+ e- We note that the only sig- 
nificant non-zero signal is contained in the {£, m) = (2, 0) 
wave mode* and we use ( [34| to get Dh+^e- When compar- 
ing the waveform obtained from CCE to the one obtained 
from RWZM (not shown), we notice that the waveforms 
from RWZM are more susceptible to numerical noise and 
contain spurious high-frequency oscillations. This is con- 



sistent with our findings in [21] (see also Sec. IIIBl. The 
waveforms extracted via RWZM are similar to those ob- 
tained in [51 which also use RWZM extraction. We 
thus believe that the results of fS', ^95] also suffer from the 
same spurious high-frequency noise. 

We align all quantities at the coordinate time when 
an apparent horizon appears (t — Ibh — 0). By com- 
puting the differences in low and medium, and medium 
and high resolutions, we get an estimate for the conver- 
gence of our simulations. In panels below the emitted 
GW signal g, and central density evolution of 



^ Earlier studies |5] |92| also found an (£, m) = (4, 0) wave mode. 
In our case, this mode is three orders of magnitudes smaller than 
the {i,m) = (2,0) mode amplitude and comparable to the level 
of numerical noise. Since the earlier study did not use causally 
disconnected outer boundaries, did not compute the waveform 
at future null infinity , and had less resolution in the wave- 
extraction zone, we argue that a (£, m) = (4, 0) could have been 
excited because of numerical artifacts and systematic errors. 



21 



TABLE V: Initial parameters and properties of the collaps- 
ing neutron star. ADM mass Madm and angular momentum 
Jadm are computed from the initial data solver at spatial 
infinity i". The radiated energy -Erad and angular momen- 
tum Jrad are computed from waves extracted via the method 
of CCE including modes up to ^ = 6. The apparent hori- 
zon mass Mah and angular momentum Jadm are computed 
on the apparent horizon surface after the black hole has set- 
tled to an approximate Kerr state. The data are reported for 
high resolution simulation r2. The value in brackets denotes 
the numerical error in the last reported digit. Units are in 
c = G = Mq = 1. 



Initial polytropic scale 




100 


Evolved polytropic scale 


K 


98 


Polytropic index 


r 


2 


Central rest-mass density 




3.116 X 10*^ 


Axes ratio 




0.65 


Rotational / binding energy [%] 


T/\W\ 


7.68 


Equatorial radius [Mq] 


Re 


9.6522 


Baryonic mass [Mq] 


Mb 


2.0443 


ADM mass [Mq] 


Madm 


1.8605 


ADM ang. mom. [Mq] 


Jadm 


1.8814 


Spin 


a 


0.5435 


Rad. energy [Mq] 


-Erad 


8.14(3) X 10" 


Rad. ang. mom. [Mq] 


J rad 


0(1) X 10^^" 


AH mass [Mq] 


Mah 


1.8602(3) 


AH ang. mom. [M|] 


Jah 


1.874(7) 



Fig. [14] respectively, we show the differences in GW sig- 
nal and central density using the three different resolu- 
tions. The differences between medium and high resolu- 
tions are scaled for second-order convergence. At black 
hole formation, the GW signal and central density exhibit 
clear second-order convergence. During collapse, while 
the central density shows second-order convergence, the 
convergence of the GW signal is somewhat obscured due 
to the oscillatory nature of the latter, especially when 
the signal is not perfectly in phase. In the lower panel of 
Fig. [14] we show the -L2-norms of the Hamiltonian con- 
straint ll-ff II2 for the three resolutions. Since the artificial 
initial pressure depletion is not constraint satisfying, the 
constraints do not converge initially. For this reason, we 
do not introduce any rescaling for convergence. However, 
the slopes for higher resolutions are smaller, resulting in 
somewhat smaller constraint violations at later times. At 
the time when an apparent horizon appears, and during 
ring-down, the constraints exhibit second-order conver- 
gence. 

In Fig. |15| we compare performance of cell-centered 
AMR with ePPM, vertex-centered AMR with oPPM, 
and cell-centered AMR with ePPM and multirate RK 
time integration using baseline resolution rl. The 
vertex-centered case with oPPM exhibits slightly larger 
constraint violations than the cell-centered setup using 
ePPM. Before the horizon forms, baryonic mass should 
be exactly conserved. In practice, this is not the case, 
even in the cell-centered case with refluxing. One reason 
for non-conservation is the artificial low-density atmo- 



sphere (see Appendix [C]) . Another reason is the buffer- 
zone prolongation in regions that thread the surface of 
the star. Here, prolongation involving cells in the at- 
mosphere can amplify mass non-conservation. We note, 
however, that the cell-centered case with refluxing per- 
forms better than the vertex-centered case. The simula- 
tion using multirate time integration performs equally 
well compared to the same simulation using standard 
RK4 time integration. 

In Fig. [16] we show the mass and spin evolution of 
the apparent horizon. After t — isH = 0, horizon mass 
and spin are quickly growing until they asymptote to- 
wards the ADM mass and angular momentum of the 
spacetime, respectively. For a given spacetime, ADM 
mass and angular momentum are always constant. Both 
quantities are calculated in the initial data solver and 
evaluated at spatial infinity. Since all matter falls into 
the horizon, the black hole mass plus the radiated en- 
ergy must be equal to the ADM mass. The same ap- 
plies to the angular momentum. In the present case, we 
have Madm = 1.8605 -M©. The black hole settles to a 
horizon mass of A/ah = 1.8602 M©. Thus, the differ- 
ence is 0.016%. Similarly, the angular momentum ini- 
tially is Jadm — 1-8814 Mq, and the black hole settles 
to J AH = 1.874 M|. This makes a difference of 0.39%. 
The radiated energy is -Erad — 8.14 x 1O~^M0 and hence 
is tiny compared to the rest mass of the system. This 
value agrees to the estimate given in jJlinil. Since the 
only significant non-zero GW mode is the (£, m) = (2, 0) 
mode, no angular momentum is radiated. We find that 
by decreasing the atmosphere level and increasing the 
resolution, the differences in horizon mass and angular 
momentum compared to the initial ADM values are de- 
creased. Hence, the error in mass and angular momen- 
tum conservation is due to systematic (atmosphere) and 
numerical error. 

In Fig. |17| we investigate the power spectrum of the 
emitted GW signal Dh+^e- The blue straight curve is 
the power spectrum of the entire signal which peaks at 
/peak = 5.06 kHz. The green dashed curve is produced 
by first applying a time-domain window function around 
the black hole ring-down part of the waveform before tak- 
ing the Fourier transform. Thus, the green dashed curve 
is the power spectrum of the black hole ring-down part 
of the waveform. This curve peaks at /peak, ring-down — 
6.47 kHz. We can compare this frequency with the the- 
oretically obtained quasi-normal (QNM) ring-down fre- 
quency for a perturbed black hole in vacuum. For the 
black hole mass A/ah = 1.8602 Af© and dimensionless 
spin a = Jah/A/^h = 0.5414, the {£,m) = (2,0) pro- 
grade fundamental (.A'^ = 0) quasi-normal frequency is 
/qmn = 6.68 kHz |129| . Thus, the relative difference 
is ~ 3.3%. This is consistent with [5] who find "good 
agreement" (unfortunately they do not provide numbers) . 
Note that we do not expect the two values to exactly co- 
incide. The theoretical QNM frequency is strictly only 
valid for perturbed Kerr black holes in vacuum. Since 
matter is crossing the horizon initially, the ring-down sig- 
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FIG. 18: Binary neutron stars: Convergence study of the 
{£, m) = (2, 2) mode of the GW strain Dh, and the L2-norm 
of the Hamiltonian constraint ||-ff||2- The top panel shows 
the "+" polarization of the (I, m) — (2, 2) mode for all three 
resolutions. The panel below shows the GW phase </> of the 
{£, m) = (2, 2) mode. The third panel from the top shows the 
difference in phase </!>, scaled for second-order convergence. 
The vertical dashed line indicates appearance of an apparent 
horizon in the high-resolution simulation. The bottom panel 
shows the L2-norm of the Hamiltonian constraint scaled for 
first-order convergence. The simulations were performed us- 
ing cell-centered AMR, refluxing, and ePPM reconstruction. 



nal will naturally be aflFected by black hole growth and 
spin-up. 



D. Binary Neutron Stars 

We investigate accuracy and convergence of the in- 
spiral and coalescence of a binary neutron star (BNS) 
system. Previous studies in full general relativity were 
restricted by the employed purely Cartesian grids (e.g. 
[6HTl] I130fll32] , also see [14| for a recent review) , and 
thus the accuracy of the GW extraction was limited. 

For the first time in the context of binary neutron star 
mergers, we use CC E for GW extraction at future null 
infinity J^^ (see Sec. II F I. This removes finite radius and 



gauge errors and, combined with our multipatch grid, 
allows us to extract the higher than leading order modes. 

Finally, we also compare vertex centered AMR with 
oPPM with cell-centered AMR with refiuxing and ePPM. 




FIG. 19: Binary neutron stars: comparison between cell- 
centered (cc) AMR with ePPM and vertex-centered (vc) AMR 
with oPPM. The top panel shows the {£, m) — (2, 2) mode of 
the "+" polarization of the GW strain Dh. The center panel 
shows the L2-norm of the Hamiltonian constraint 11-^112. The 
bottom panel shows conservation of baryonic mass AIb ■ The 
vertical dashed line indicates the appearance of an apparent 
horizon in the baseline resolution simulation. The simulations 
were performed using resolution rl, though for the conserva- 
tion of mass, we also show the high resolution (r2) result. 
The error in mass conservation converges with better than 
second-order as the resolution is increased up to the point 
when a new refinement level is switched on at t ~ 7.5 ms. 



1. Initial Conditions and Equation of State 

The particular system we evolve is the initial data 
set G2_I12vsl2_D5R33_60km produced by the LORENE 
code |57[I133'. This system, with the same parameters as 
described below, has also been considered in [134, ,il35] . 

The system consists of two neutron stars initially de- 
scribed by a polytropic equation of state P = Kp^ with 
K — 123.6 and F = 2 with an initial coordinate sep- 
aration of 45 km. We evolve the system using a F-law 
equation of state of the form 



P=(r-l)pe. 



(35) 



These parameters yield neutron stars of individual bary- 
onic mass Mb — 1.78 M© and ADM-mass in isolation 
of A/ns = 1.57 Mq. The total ADM mass of the sys- 
tem is A/adm = 3.2515 Af0, and the total ADM angular 
momentum is Jadm — 10.1315 Afg. The initial orbital 
angular frequency of the binary is iljni = 302 Hz. The 
initial parameters and properties are listed in Table |VI[ 



2. Numerical Setup 

The numerical setup consists of the six spherical 
inflated- cube grids that surround the central Cartesian 
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FIG. 20: Binary neutron stars: GW modes {£, m) = 
(3,2), (4,4), (6,6), (8,8) of "+" polarization of the strain Dh 
unambiguously extracted via CCE. The waveforms are shown 
for high resolution simulation r2. The vertical line indicates 
the time of appearance of an apparent horizon. Following 
appearance of an apparent horizon, a black hole ring-down 
signal is visible. 



cube. The inner spherical radius of the inflated cube 
grids is located at a coordinate radius of i?s — 75.84 M0 
and the outer (spherical) boundary is located at 
dius of i?B — 2800 A/0. The radial resolution at the 
inner spherical inter-patch boundary matches the coarse- 
grid Cartesian resolution of the central cube and is 
Ax = 1.5 Mq = 2.22 km, Ace = 1.2 M© = 1.77 km and 
Aa;0.96 Mq = 1.42 km for the low, medium and high res- 
olution runs, respectively. In the region 250 < r < 
800 M0 we smoothly transition to a coarser resolution 
of 6.OM0, 4.8M0 and 3.84M0 for low (rO), medium 
(rl) and high resolution (r2), respectively. The angular 
resolution is constant along radial distances and we use 
21, 25 and 31 angular grid points per angular direction 
and spherical patch for the three resolutions. We use 
4 initial levels of mesh reflnement in the inner Cartesian 
cube to resolve the neutron stars. We surround each neu- 
tron star with a set of nested, refined cubes of half-width 
13 M0, 17.875 M0 and 26.125 M0, where the finest level 
completely covers the neutron star. All refined cubes sur- 
rounding the stars are contained in the common, coarse 
cube of half- width R^. In each refined level the resolu- 
tion is twice that of the previous level. On the finest 
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FIG. 21: Binary neutron stars: Phase convergence of the 
(£, m) = (6, 6) mode of the GW strain Dh. The top panel 
shows the "+" polarization component Dh^ , and the panel 
below shows the phase (j), for low rO, medium rl, and high 
r2 resolutions. The bottom panel shows the phase differ- 
ences between low and medium, and medium and high reso- 
lutions, scaled for second-order convergence. Convergence is 
maintained throughout inspiral and merger. In the ring-down 
phase, the coarsest resolution rO is insufflcient to accurately 
resolve this mode, and the results cease to converge properly. 
The vertical line indicates appearance of an apparent horizon 
in the high resolution simulation. 



level, the neutron stars are covered with a resolution of 
Aa; = 0.1875 Af0 = 0.278 km, Ax = 0.15 M0 = 0.222 km 
and Aa; = 0.12 M0 = 0.176 km for the three resolutions 
rO, rl and r2, respectively. 

When the two neutron stars are about to come into 
contact, we remove the nested set of cubes surrounding 
each individual star and surround the binary with a com- 
mon set of nested cubes of half- width i?s, 30 M0, 15Af0 
and 7.5 Af0 ensuring uniform resolution in the central 
region. Once the lapse function drops to values that in- 
dicate that an apparent horizon is about^ to form, we 
switch on a final level of radius 3.5 Af0 and resolution 
9.38 X 10-2 Af0, 7.5 x IQ-'^ Mq and 6.00 x 10"^ Af0 for 
the low, medium and high resolution runs respectively. 
This level allows us to handle the steep metric gradients 
developing inside of the newly formed apparent horizon. 

During inspiral, we track the center of mass of each 
neutron star to keep the two fiuid bodies close to the 
center of their refined regions. We compute the center 
of mass of an individual neutron star by integrating over 
the conserved density within a radius R = 4.0 Af0 of 



^ This is a consequence of the 1 -|- log slicing condition which 
locally slows down time evolution (i.e. a < 1) in regions of strong 
curvature. A closed surface of lapse of a < 0.3 has been found 
to approximately resemble the apparent horizon shape. 
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FIG. 22: Binary neutron stars: Amplitude convergence of 
the {£, m) = (6, 6) mode of the GW strain Dh = Dh{(j)) as 
a function of phase cf). The top panel shows the "+" polar- 
ization component Dh^{cf>), and the panel below shows the 
amplitude A{(j>), for low (rO), medium (rl), and high (r2) 
resolutions. The bottom panel shows the amplitude differ- 
ences between low and medium, and medium and high reso- 
lutions, scaled for second-order convergence. Convergence is 
maintained throughout inspiral and merger. In the ring-down 
phase, however, the coarsest resolution rO is insufficient to ac- 
curately resolve this mode, and the results cease to properly 
converge. The vertical line indicates appearance of an appar- 
ent horizon in the high resolution simulation. 



the densest point on the grid. This method produces 
smoother tracks than directly using the location of the 
densest point, and helps reducing the jitter in the mesh 
refinement boxes observed otherwise. 

We set the damping coefficient of the F-driver gauge 
condition to = 1. 

We set the dissipation strength to ediss = 0.1 every- 
where on the grid. The artificial low-density atmosphere 
is 10^ times lower than the initial central density. 



3. Discussion 

While the two neutron stars orbit each other, they lose 
energy due to gravitational radiation, inspiral, and finally 
merge. The nascent hypermassive neutron star remnant 
has a mass which is well above the maximum mass of 
neutron stars. It forms a black hole on a dynamical 
timescale. The black hole is initially highly excited, and 
relaxes to a Kerr state by emitting gravitational ring- 
down radiation. 

In Fig. [18] we show convergence of the dominant 
{£, m) — (2, 2) mode of the GW strain Dh, the i2-norm 
of the Hamiltonian constraint ||i/||2- The upper panel 
shows the "+" polarization of the (€, m) — (2, 2) mode of 
the GW strain for the resolutions rO. rl, and r2. The 



TABLE VI: Parameters of the binary neutron star system. 
ADM mass A/adm and angular momentum Jadm are com- 
puted by the initial data solver at spatial infinity i". The 
radiated energy -Brad and angular momentum J^ad are com- 
puted from waves extracted via CCE including modes up to 
I = Q. The apparent horizon mass Mah and angular mo- 
mentum Jadm are computed after the black hole has settled 
to an approximate Kerr state. Gravitational disk mass Afdisk 
and angular momentum Jdisk are calculated from energy and 
angular momentum conservation. The data are reported for 
simulation r2. The value in brackets denotes the numerical 
error in the last reported digit. Units are in c = G = Mq = 1. 



Lorene initial data set 

Initial separation [km] 

Polytropic scale 

Polytropic index 

Initial orbital frequency [Hz] 

ADM mass [Mq] 

ADM ang. mom. [Mq] 

Rad. energy [Mq] (%) 

Rad. ang. mom. [Mq] (%) 

AH mass [Mq] 

AH ang. mom. [Mq] 

AH spin 

Grav. mass disk [Mq] 
Bary. mass disk [Mq] 
Ang. mom. disk [Mq] 



G2 

d 
K 

r 

Madm 
Jadm 

-E-rad 
J rad 

Mah 
Jah 
a 

Mdisk 

Ms, disk 
J disk 



I12vsl2_D5R33_60km 
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123.6 
2 

302 
3.2515 
10.1315 

2.51(5) X IQ-^ (0.77%) 
1.206(9) (11.9%) 
3.2249(3) 
8.75(2) 
0.841(2) 
1.4(4) X 10^^ 
1.3(2) X 10"^ 
0.16(4) 



waveform is extracted via CCE. To obtain Dh, we use a 
cut-off parameter /q = 507 Hz, which is below the initial 
instantaneous {t, m) — (2, 2) mode frequency f?^^ deter- 
mined from the initial orbital frequency by = 20.^^1. 
To assess the phase convergence, we plot the differences 
in phase between low rO and medium rl resolution, and 
medium and high r2 resolution, scaled for second-order 
convergence. We also plot the L2-norm of the Hamilto- 
nian constraint ||i?||2 scaled for first-order convergence. 

ITTlSI 



Similar to the isolated neutron star tests in Sec 
the dominant constraint error is generated at the con- 
tact discontinuity at the neutron star surface, where our 
scheme locally reduces to first-order accuracy. 

In Fig. 19 we compare cell-centered (cc) AMR and 
ePPM reconstruction with vertex-centered (vc) AMR 
and oPPM. The {I, m) = (2, 2) mode of the GW strain 
Dh and the L2-iiorm of the Hamiltonian constraint ||i?||2 
do not show any significant differences between the two 
numerical setups at this point. After black hole and disk 
formation, the vertex-centered scheme exhibits a slightly 
larger slope in constraint growth. In the bottom panel, 
we show conservation of total baryonic mass Mb- Dur- 
ing early inspiral, both setups conserve mass to a high 
degree, only affected by small errors due to our artifi- 
cial atmosphere (see Appendix [C]) . Note that both neu- 
tron stars are completely contained on their finest grids. 
Thus, there are no refinement boundaries directly infiu- 
encing the evolution of the two fiuid bodies. As the in- 
spiral progresses, we find that mass conservation is vio- 
lated in the cell-centered case to a higher degree than in 
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the vertex-centered case (though the error converges as 
the resolution is increased). This appears to be an arti- 
fact of buffer zone prolongation close to the neutron star 
surface in combination with low density matter slightly 
above and at atmosphere values. Due to numerical er- 
rors, small amounts of mass are leaking out of the neutron 
star during inspiral and interact with the atmosphere. As 
this low density matter reaches the buffer zones, numeri- 
cal errors due to prolongation, which are by construction 
larger in the cell-centered case, tend to amplify the nega- 
tive effects of the atmosphere treatment. In experiments 
with isolated neutrons stars, however, we find that when 
the refinement boundaries are sufficiently far removed, 
and/or the atmosphere level is further decreased, mass 
can be conserved to a higher degree. 

We also compare the simulations to a setup using mul- 
tirate RK time integration and cell-centered AMR with 
ePPM. Unfortunately, due to the large fluid bulk veloc- 
ities in the inspiral phase, the orbital phase accuracy is 
significantly affected by the lower order fluid time in- 
tegration. Thus, we do not recommend application of 
multirate RK schemes in the context of binary neutron 
star mergers, especially when orbital phase accuracy is 
paramount. The problem may be ameliorated by the use 
of co-rotating coordinates (see, e.g., |71j). 

In order to demonstrate the potential of the multipatch 
scheme for more accurate wave extraction, we show in 
Fig. [20| some of the higher harmonic GW modes that are 
emitted during inspiral, merger, and ring-down. We show 
(from top to bottom) the {£,m) ^ (3,2), (£,m) = (4,4), 
(£,m) = (6,6), and {i,m) = (8,8) modes of "+" polar- 
ization of the strain Dh. The modes are extracted from 
a simulation using resolution r2, cell-centered AMR, and 
ePPM. All modes up to {£, m) = (4, 4) show a clean in- 
spiral, merger and ring-down signal, and converge with 
resolution (see below). For higher modes, our lowest 
resolution run rO is insufficient to also allow for clean 
convergence of the corresponding ring-down signals. Ac- 
cordingly, those should be taken with a grain of salt. As 



an example, in Figs. 21 and |22[ we show convergence of 
phase and amplitude of the (£, to) — (6, 6) mode of the 
GW strain, respectively. Fig. |2T] shows the GW ampli- 
tude A reparametrized in terms of the gravitational phase 
(f) to disentangle phase from amplitude. Both figures indi- 
cate that second-order convergence is maintained during 
inspiral up to merger. The ring-down part, however, does 
not exhibit clean second-order convergence. In that case, 
the coarse resolution becomes insufficient, and the result 
ceases to converges properly. We note that for the highest 
extracted mode, {£,m) — (8,8), the coarsest resolution 
is insufficient to allow for clean convergence also in the 
inspiral phase. 

We compute the radiated energy i^radj radiated angu- 
lar momentum Jiad, the horizon mass Mahi and horizon 
angular momentum J ah- For the computation of the ra- 
diated quantities, we include modes ^ ^ 6 as extracted 
via CCE. After the black hole has formed and settled to 
an approximate Kerr state, some amount of material is 



located in an accretion disk surrounding the black hole. 
Hence, we do not expect that horizon mass and radiated 
energy balance with the total ADM mass at this time. 
Rather, the difference denotes the gravitational mass of 
the accretion disk that has formed. Likewise, the same 
is true for the balance of angular momentum. Given the 
horizon mass, the spacetime's total ADM mass, and the 
radiated energy, we estimate the gravitational mass of 
the accretion disk to be Mdisk = -^^adm — Mah — ^-rad = 
(1.4 ± 0.4) X 10~^ Mq. The disk's baryonic mass is 
A^s.disk = (1-3 ± 0.2) X 10^^ Mq, which we compute by 
integrating over all material outside of the apparent hori- 
zon and within a radius R < 40 M0. Both, baryonic and 
gravitational mass agree within their error bars. We note 
that the mass of the disk, though clearly visible in den- 
sity contour plots of our simulation (not shown), is tiny 
and thus not much above the numerical error. Given the 
horizon angular momentum, the spacetime's total ADM 
angular momentum, and the radiated angular momen- 
tum, we estimate the disk's angular momentum to be 



0.16 ± 0.04 M|. For conve- 
nience, we list spacetime, black hole, disk, and radiated 



i^disk — J A 



DM - 



Jah — Ji 



rad 



mass (and angular momentum) in Table VI All error 



bars are estimated using medium and high resolution re- 
sults. The results for mass and spin of the black hole 
agree to the values that were found in |134| . 

In our binary neutron star merger problem, we also in- 
vestigate the error inherent to finite-radius GW extrac- 
tion. We compare ^'4 as extracted via the NP formalism 
at a finite radius with \l/4 as extracted via CCE at future 
null infinity J'^ . We align two given waveforms in the 
early inspiral phase by minimizing their phase difference 
over an interval t Cz [2.5 ms, 3.5 nis] using the method de- 
scribed in [136 . For the {£,m) — (2,2) mode, we find 
a total dephasing on the order of A(f) ^ 1 rad and an 
amplitude difference of about ~ 10% between the wave- 
form obtained at i? = 25OM0 and the one obtain at 
J^'^ . Waveforms extracted at smaller radii naturally yield 
larger differences to the result at J'^ . While the ampli- 
tude error is rather large, the dephasing is comparable 
to the dephasing due to numerical error of the orbital 
evolution of the two neutron stars. Since this numerical 
error is convergent, but the systematic finite radius-error 
is not, the finite-radius error becomes a non-neglegible 
effect as the numerical resolution is increased. As shown 
in ^T, '23] for the case of binary black hole mergers, ex- 
trapolation to infinity using finite-radius data can reduce 
the errors to a tolerable level in cases where CCE is not 
available. 

Finally, we investigate the influence of the outer 
boundary when it is not causally disconnected from the 
wave extraction region and interior evolution. We com- 
pare a setup with a causally connected outer boundary 
located at i?B = 2OOOAf0 and a causally disconnected 
boundary located at Rb — 28OOM0 . The former setup is 
in causal contact with the interior and wave-extraction 
region during the merger and ring-down phases. We find 
a difference in GW phase and amplitude, and final spin 
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and mass of about ~ 7%. More details are given in Ap- 
pendix ^ 

By comparing our results with those of |134[ 1135] , we 
conclude that the accuracy of the orbital evolution of the 
two neutron stars is very similar. The errors in satisfying 
the Hamiltonian constraint and conserving baryonic mass 
are of comparable size. This is not surprising, since we 
find little difference between the new cell-centered AMR 
scheme compared with the vertex-centered AMR scheme 
that was also used in |134l 1135] . Due to our multipatch 
grids, causally disconnected outer boundaries, and CCE, 
however, the waveforms that are extracted from our sim- 
ulations are more accurate than what has been shown in 
previous studies. 



IV. SUMMARY AND CONCLUSIONS 

We have presented a new GR hydrodynamics scheme 
using multiple Cartesian/curvi-linear grid patches and 
flux-conservative cell-centered adaptive mesh refinement 
(AMR) to allow for a more efficient and accurate spatial 
discretization of the computational domain. This is the 
first study enabling GR hydrodynamic simulations with 
multipatches and AMR. Our multipatch scheme consists 
of a set of curvi-linear spherical "infiated-cube" grids with 
fixed angular resolution and variable radial spacing, and 
a central Cartesian grid with AMR. High-order Lagrange 
interpolation is used to fill ghost zones at patch bound- 
aries for variables that are smooth, and second-order 
essentially non-oscillatory (ENO) interpolation for vari- 
ables that contain discontinuities and shocks. 

Apart from the successful implementation of mul- 
tipatches and flux-conservative cell-centered AMR, we 
have introduced a number of additional improvements 
to the publicly available code GRHydro: (i) We have ap- 
plied the enhanced piecewise-parabolic method (PPM) 
to ensure high-order reconstruction at smooth maxima, a 
property that we have found to be crucial for cell-centered 
AMR. (ii) To speed up the computation, we have applied 
a multirate Runge-Kutta time integrator that exploits 
the less restrictive Courant-Friedrich-Lewy (CFL) condi- 
tion for the hydrodynamic evolution by switching the the 
time integration to second order and thus reducing the 
number of intermediate steps by a factor of two. Since the 
hydrodynamic evolution dominates the curvature evolu- 
tion in terms of computational walltime when complex 
microphysics and neutrinos are included, the scheme can 
yield a speedup of > 30% (e.g. [90i). 

We have presented stable and convergent evolutions 
for binary neutron star mergers, stellar collapse to a neu- 
tron star, neutron star collapse to a black hole, and evo- 
lutions of isolated unperturbed and perturbed neutron 
stars. For each test case, due to the more efficient domain 
discretization, we have been able to enlarge the domain 
sufficiently so that the outer boundary is causally discon- 
nected from the interior evolution and wave-extraction 
zone. This has allowed us to remove the systematic error 



that arises from the lack of constraint preserving bound- 
ary conditions for the Einstein equations in the BSSN for- 
mulation. In the case of the binary neutron star merger 
problem, we have found that this error is on the order of 
a few percent, and thus limits the accuracy of the simu- 
lation and GW extraction. 

In addition to enlarging the domain, multipatches have 
also allowed us to significantly increase the resolution in 
the GW extraction zone compared to previous studies. 
For the neutron star merger problem, we have been able 
to extract convergent spherical harmonic modes of the 
GW strain Dh up to ^ = 6. Previous studies have only 
considered the dominant {£, m) — (2, 2) wave mode for 
this problem. 

Furthermore, we have been able to remove the sys- 
tematic error inherent in finite-radius wave extraction by 
application of Cauchy-characteristic extraction (CCE). 
This wave-extraction method computes gauge-invariant 
radiation at future null infinity ^7+ using boundary data 
from a worldtube at finite radius. This method has pre- 
viously been applied in simulations of binary black holes 
and stellar collapse ^2^ ,21n24l l75l [76] . Here, we have 
applied CCE also to simulations of binary neutron star 
mergers, neutron star collapse to a black hole, and iso- 
lated excited neutron stars. We have found that the error 
due to finite-radius extraction can be as large as 10%. 

Finally, for each test case, we have compared the orig- 
inal vertex-centered AMR scheme using original PPM 
with the new fiux-conservative cell-centered AMR scheme 
using enhanced PPM. The accuracy has been investi- 
gated and compared to results from previous studies. We 
have found that simulations of stellar collapse greatly 
benefit from fiux-conservative cell-centered AMR with 
enhanced PPM compared to the original vertex-centered 
AMR scheme with original PPM. Conservation of mass 
and the satisfaction of the Hamiltonian constraint are 
significantly better with the new scheme. The isolated 
neutron star and binary neutron star test cases, on the 
other hand, are not much affected by the choice of cell- 
centered or vertex-centered AMR. This is mainly due to 
the choice of grid setup: no matter is crossing any re- 
finement boundaries so that fiux-conservation is not im- 
portant. It can become important, however, in the post- 
merger phase of binary neutron star coalescence, espe- 
cially in cases where a massive accretion torus forms. 

The multipatch infrastructure, the associated curva- 
ture and hydrodynamics evolution codes, and all other 
computer codes used in this paper will be made (or are 
already) publicly available via the EinsteinToolkit [91j. 
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FIG. 24: Conservation of mass (top panel), energy (middle 
panel), and momentum (bottom panel) as a function of time 
for a shock front crossing a refinement boundary. The solid 
(red/blue/green) lines are from a simulation with refluxing, 
while the dashed (black) curves show the case without re- 
fluxing. With refluxing, mass, energy, and momentum are 
exactly conserved (to machine precision). Without refluxing, 
conservation of mass, energy, and momentum is violated. 



Appendix A: Shock-tube Tests 

We perform a number of basic Sod shock tube and 
spherical blast wave tests on fixed backgrounds to en- 
sure correctness and convergence of our scheme at mesh- 
refinement and inter-patch boundaries. 

In this appendix, we restrict our attention to a simple 
Sod test to show convergence of the primitive variables 
across inter-patch boundaries (see Sec. IIC3I, and to 



demonstrate mass, energy, and momentum conservation 
at refinement boundaries when refiuxing (see Sec. II D I 
is used. 

The Sod shock-tube test consists of setting the initial 
fluid state according to |137| . The shock front is located 
at a position xq. The background metric is set to the flat 
space Minkowksi metric. The tests below use a gamma- 
law equation of state P = {T — l)pe with T = 1.4. 

If not stated otherwise, the tests below use cell- 
centered AMR with refluxing, ePPM reconstruction, 
second-order ENO inter-patch interpolation, RK4 time 
integration with At/ Ax — 0.4, and the HLLE Riemann 
solver. 



1. Inter-patch Interpolation 

In this particular test, we check that shock fronts are 
correctly transported across inter-patch boundaries by 
maintaining convergence, and without introducing local 
oscillations at the shock, even in the presence of non- 
trivial Jacobians and coordinate transformations. We 
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setup a multipatch grid consisting of a central Cartesian 
grid surrounded by the spherical inflated-cube grids. The 
outer boundary extends to Rb = 2.5Mq. The bound- 
ary between Cartesian and spherical grids is located at 
Rs = O.5M0. No AMR is employed. For the coarsest 
resolution (rO), we set the Cartesian and radial resolu- 
tion to Ax = Ar = 0.05, and use {Np,N„) (20,20) 
cells per spherical patch per direction. Medium (rl) and 
high (r2) resolutions double and quadruple, respectively, 
the resolution with respect to the coarsest resolution. 

We set Sod initial data with a;o = and evolve the 
system for sufficiently long so that the shock propagates 
across inter-patch boundaries. At each timestep, we com- 
pare the evolved fluid state with a solution from an exact 
special relativistic Riemann solver |138| . 

In Fig. [23] we show the ii-norm of the difference be- 
tween exact and evolved primitive density p, speciflc in- 
ternal energy e, and the x-component of the 3-velocity 
v^. All quantities are plotted for the three resolutions 
rO, rl, and r2. As the resolution is increased, the error 
correctly decreases by a factor of two between successive 
resolutions, thus indicating first-order convergence. This 
is consistent with the ENO operator, which reduces to 
first-order at shocks. 
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III A[ ): the effect of original 
ePPM) on the Hamil- 
tonian constraint as a function of x at time t — 0.76 ms. on 
cell-centered (cc) and vertex-centered (vc) AMR grids. The 
star's radius is Re = 14.16 km. The original PPM results 
in large constraint violations on the cell-centered grid. The 
enhanced PPM clearly outperforms oPPM. For ePPM, the 
error is dominated by the first-order error at the neutron star 
surface, where the scheme reduces to first order. 



Appendix B: Enhanced PPM Scheme 



2. Refluxing 

In this simple test, we check the correctness of our re- 
fiuxing scheme with a shock front crossing a refinement 
boundary. As the shock crosses the boundary, mass, mo- 
mentum and energy must be conserved to machine pre- 
cision. 

The numerical grid consists of two levels of 2:1 AMR. 
The coarse level extends from x ~ to x = 1. The fine 
level has a refinement half- width of r — 0.1 and is located 
at X = 0.4. We set the Sod shock front |137| at location 
xq = 0.48. Thus, the shock starts off on the fine grid and 
propagates onto the coarse grid. 

A measure of conservation of energy and mass is given 
by the sum of the conserved internal energy r and the 
conserved density D over the entire simulation domain, 
respectively. Both sums must be constant for all times 
t. A measure for conservation of momentum is given by 
the balance between the conserved momentum and the 
pressure force per unit time. The balance as a sum over 
the entire simulation domain must be constant as a func- 
tion of time. In Fig. [24] we show the sums of conserved 
density, energy, and momentum when refiuxing is used 
(solid lines). Without refluxing (dashed lines), the con- 
served mass, energy, and momentum grow significantly 
at time t ~ 0.025 when the shock front crosses the refine- 
ment boundary. 



The PPM scheme seeks to find "left" and "right" in- 
terpolated values, ai^L and Ui^R at the left and right cell 
interfaces of a primitive quantity defined on cell cen- 
ters labeled by i = 0, ..,N — 1. The left and right states 
are defined on cell interfaces labeled by a^j.i. Rather 
than assuming a constant value for a cell-averaged quan- 
tity within a given cell, the PPM scheme uses parabolas 
to represent cell averages within a given cell. 

The enhanced PPM reconstruction proceeds in three 
steps: (i) Compute an approximation to a at cell in- 
terfaces using a high-order interpolation polynomial, (ii) 
limit the interpolated cell-interface values obtained in (i) 
to avoid oscillations near shocks and other discontinu- 
ities, (iii) constrain the parabolic profile so that no new 
artificial maximum is created within one single cell. The 
main difference to the original PPM scheme is in steps 
(i) and (ii). Both, the limiter and the constraining of 
the parabolic profiles is more restrictive in the original 
PPM scheme, thus reducing the order of accuracy in cases 
where it is not necessary. 

a. First Step: Interpolation We compute an approx- 
imation to a at cell interfaces, which, assuming a uniform 
grid, is obtained via fourth order polynomial interpola- 
tion 

7 1 

using the cell center values of a from neighboring cells. 
Ref. |86| also suggests to use a sixth-order polynomial. 
This, however, requires more ghost points. In our tests, 
we find no significant difference between fourth and sixth- 
order interpolation. Hence, we stick to the fourth-order 
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interpolant. 

b. Second Step: Limiting We require that the values 
aj_|_i satisfy 



(B2) 



i.e., the interpolated value a^_^_l must lie between adja- 
cent cell values [86 . This is enforced by the following 
conditions. If (B2| is not satisfied, then we define the 



second derivatives, 



3(aj - 2a^^i + a,+i) 



(B3) 



{D^a)^^i^ := (a,_i - 2aj + a,+i) , (B4) 



(a^ - 2ai+i + ai+2) ■ 



(B5) 



If (I?^a)j_|_i and {D^a)^_^_l ^ all have the same sign s = 
sign((Z3^a)j^i ), we further define 

(.D'^a),+ i^ii„, := smm{C\{D'^a)^^i j^\, 
C\iD^a),+ .,j,\, 
\{D'a)^^.\). (B6) 

where C is a constant that we set according to [86J to 
C — 1.25. Otherwise, if one of the signs is different^, we 
set {D'^a)i^i lijjj = 0. Then, we recompute (Bll by 



1 



1 



fli+i = ^{a^ + fli+i) - -{D a),,+ i 



,lim ' 



(B7) 



c. Third Step: Constrain Parabolic Profiles Here, 
we apply the refined procedure from [87] . We begin by 
initializing left and right states according to the interpo- 
lated (and possibly limited) a^_^i via 



ai,R — ai+i,L 



(B8) 



so that the Riemann problem is trivial initially. The 
conditions below potentially alter 7^ and a^+i i, so that 
the Riemann problem becomes non-trivial. 

First, we check whether we are at a smooth local max- 
imum. A condition for local smooth maxima is given by 



(ai,L - ai)(ai - at^R) < , 
{ai-2 - ai){ai - ai+2) < 0. 



or 



(B9) 



For the specific internal energy e, we also set {D^a)^_^_i = 0, 

in cases when {D^a)^_^_i > ^(fii -|- ai_|_i). This is different 

from the procedure in 1871 , but is necessary at very strong contact 
discontinuities such as the surface of a neutron star to prevent e 
from becoming negative for equations of state that do not allow 
e < 0. In practice, this additional limiter has no effect on the 
measured accuracy. 



If (B9l holds, we compute, similar to (B3l, 

{D'^a)i^L = 0,1-2 - 2aj_i -t- at , 
[D'^a)i^R ^ at - 2ai+i + ai+2 ■ 



(BIO) 



If (-D^a)j jj] all have the same sign s — sign{{D^ a) i) , 
we compute 

(£'^a),+ i = smm{C\iD^a),+ i^L\, 

|(Z?H+il)- (Bll) 

Otherwise, if one of the signs is different, we set 
(Z?2a)^^, j.^ = 0. If 

\{D'^a)i\ < 10"^^ • max(|ai_2|, |ai-i|, \ai\, \ai+i\, \ai+2\) 

(B12) 



then we define and set pi — 0. Otherwise, we define 



Pi 



(B13) 



To avoid limiting at small oscillations induced by round- 
off errors, we do not apply any limiter if > 1 — 10"^^. 
Otherwise, we compute the third derivative according to 



We set 



3 _\min 



and 



3 _\max 



= (^'a)^+l,c - (i?'a)^,c ■ 



max((D3a),^_3,(D3a)^_i, 



Then, we test if 

C3-max(|(i?3^)r'^|,|(i?3„)--|) 



3 _\max 



3 \iniii 



(B14) 



(B15) 



(B16) 



(BIT) 



holds. In the expression above, C3 = 0.1, according to 



Ref. [87 . If (BIT I does not hold, a limiter is not applied. 
Otherwise, we test the following conditions: (i) if (a^^^ — 
a.i){ai — Qi^fi) < 0, we set 

ai,L = ai- pi{ai - Cj^l) , 

ai,R = ai + pi{ai^R - Oi) . (B18) 

Otherwise, (ii) if [a^ — a^.^l > 2|ai,j^ — a^l, we set 

ai,L = tti- 2(1 - pi){ai^R - at) - pi{ai - at^L) (B19) 
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or (iii) if |oi,_r - at] > 2\ai - ai^L\, we set 

fli.fl = fli + 2(1 - pi){ai - a.ix) + Pi{ai,B^ ~ cii) . (B20) 

In the conditions (i)-(iii) above, we introduce a spe- 
cial treatment for the specific internal energy e. If 
\ai - Ui^il < |oi| or la^^R - ai\ < {a^l, we set ai,L,_R = 
instead of using the full expressions, respectively. This is 
different from the original procedure of Ref . [ST] . It es- 
sentially reduces the reconstruction of e to first order in 
cases when the correction becomes larger than the value 
of the reconstructed quantity itself. This is similar to the 
limiter step further above and is necessary at very strong 
contact discontinuities such as the surface of a neutron 
star. Without this additional limiter, e may become ill- 
conditioned. This typically happens when e is very small 
and the correction becomes larger than e itself poten- 
tially leading to negative e. For some equations of state, 
e < is ill-defined, causing the HLLE Riemann solver to 
fail. In practice, this reduction does not affect the overall 
accuracy of the scheme. We also note that this special 
treatment does not forbid e from becoming negative. 
Finally, we recompute a^.i according to 



0'i,L{R) 



a-i) 



lim 



(B21) 



In case the denominator becomes zero in the expression 
above, we set the last term to zero. 

Finally, if (B9l does not hold, we test whether 
Wi.R{L) ~ > '2\o-i,L(R) ~ o.i\ holds. In that case, we 
set 



,R(L) = flj - 2(a,j^ 



L{R) 



a-i) 



(B22) 



for either o^^l or ai^R, respectively. In the case of recon- 
structing the specific internal energy e, if |ai — 2aix{R.) I > 
Oil, we simply set a^ jK^i^^ = ai. This is for the same rea- 
son that has been mentioned above already. 

After having obtained a^^^ and ai^R, we apply the 
"standard" fiattening procedure discussed in the Ap- 
pendix of [97 . This completes the enhanced PPM scheme 
applied in our code. Note that Ref. [57] (in contrast to 
[86J) suggests to skip the second step. In our experiments 
with an excited neutron star and a collapsing stellar core, 
however, we find that when skipping this step, the scheme 
becomes too dissipative. 

The enhanced PPM scheme requires four ghost points. 
For efficiency reasons, it may be desirable to use only 
three ghost points, since less memory and interprocessor 
communication is required. In order to reduce the num- 
ber of required stencil points to three, we use fourth-order 
polynomial interpolation (Bll instead of sixth-order in- 
terpolation [SB] in the first step, and we skip the check 
(B17) involving the third derivatives {D^a)i. We also use 



a modified fiattening scheme which allows us to use only 
three ghost points. This modified fiattening scheme is 
the same as the one presented in the Appendix of ^97j, 



but we drop the maximum in Eqn. (A. 2) of [HZj, and di- 
rectly use fi = fj. In our tests, we have found only small 
differences between the four- and three-point scheme. 

In Fig. [25] we show the effect of ePPM compared with 
oPPM on the Hamiltonian constraint H along the x-axis 



for the example of an isolated TOV star (Sec. Ill A I on 
cell-centered and vertex-centered AMR grids. Clearly, 
ePPM results in a significantly lower error compared to 
oPPM on vertex-centered, and especially on cell-centered 
AMR grids. 



Appendix C: Atmosphere Treatment 

In vacuum, obviously, the equations describing the 
fiuid dynamics break down. When simulating isolated 
neutron stars or binary neutron star mergers, a large frac- 
tion of the simulation domain is physically vacuum. At 
the surface of the fiuid bodies where a sharp transition 
to vacuum occurs, the Riemann solver breaks down. 

As a simple solution to this problem, we keep a very 
low and constant density fiuid (the atmosphere) in the 
cells which would be vacuum otherwise. We also keep 
track of where the evolution of the fiuid variables fails to 
produce a physical state and reset these cells to atmo- 
sphere. Typically, there are few such cells, which cluster 
around the surface of the star. The atmosphere density 
Patmo is usually chosen to be 8 to 10 orders of magnitudes 
lower than the central density of the fiuid body. This en- 
sures that the atmosphere does not contribute noticeably 
to the total rest mass and energy in the simulation. 

Whether a given fiuid cell is set to atmosphere val- 
ues is decided depending on the local fiuid density. If it 
drops below atmosphere density Patmo, the cell is set to 
atmosphere density with zero fiuid velocity. 

More specifically, we proceed in the following way. 

1. During each intermediate time step, we set an 
"atmosphere" fiag in an atmosphere mask Ma if 
T + AtRr < or D + AtRo < 0, where i?r and Rd 
are the right-hand sides of the r and D equations 
([5|, respectively and At is the temporal timestep 
size. In addition to setting the atmosphere fiag, 
we also set all fiuid right-hand sides for that cell 
to zero, in effect freezing the further evolution of 
this cell. In that case, we also skip conversion of 
conserved to primitive variables of that cell. 

2. After a full time step, we set all variables of those 
cells to atmosphere values that are flagged as at- 
mosphere. 

3. Finally, we clear the atmosphere mask M^. 

Furthermore, we perform the following operations in- 
volving atmosphere checks: 

1. After reconstruction, we check whether the recon- 
structed primitive density is below atmosphere den- 
sity. If this is the case, we enforce first order 
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reconstruction, i.e. we set left and right cell face 
<ii,L = o,i,R = to the cell average for all prim- 
itive variables. 

2. At the end of conservative to primitive conversion, 
we check whether the new set of primitive variables 
is below atmosphere level for a given cell. If this is 
the case, we reset that cell to atmosphere level. 

In the two cases above, the atmosphere mask is not set. 

To limit high-frequency noise in cells slightly above 
atmosphere level, we set cells to atmosphere value if they 
are within a given tolerance 5 above atmosphere density, 
i.e. we test whether 



< 



Patmo (1 + (5) • 



(CI) 



In the cases considered here, we set 5 = 0.001. 

The particular treatment of vacuum regions by enforc- 
ing a low density atmosphere is not ideal and has several 
drawbacks. If a cell is forced to be not lower than a 
particular minimum density, small amounts of baryonic 
mass can be created or removed. This breaks the strictly 
conservative nature of our hydrodynamics scheme and 
can thus lead to small errors. As noted in |139| . intro- 
ducing an artificial atmosphere may also change the local 
wave structure of the solution. An artificial low density 
atmosphere can be avoided by modifying the Riemann 
solver at those cells adjacent to vacuum cells |il39j . In 
practice, however, if the atmosphere level is sufficiently 
low, the negative influence on the fluid evolution can be 
neglected. 



Appendix D: Scheduling of Ghost-Zone 
Synchronization 

We find that excessive inter-processor and inter-patch 
synchronization of ghost zone information can lead to 
significant performance drawbacks, especially on large 
numbers of processing units (> 1000). We have thus op- 
timized our ghost-zone update pattern and reduced the 
number of necessary synchronization calls. 

We distinguish between three different synchronization 
update operations: (i) inter-processor and inter-patch 
synchronizations performed after each intermediate time 
step, and (ii) AMR buffer-zone prolongation performed 
after each full time step, and (iii) AMR prolongation after 
regridding (see [81 on the latter two cases for details). 

We distinguish between two sets of vari- 
ables. One set is comprised of the spacetime 

variables 1 0, , K, Aij , f % a, /3* , | 

the curvature evolution and gauge 
and the other set is comprised 
{D, T, S,,p, e, v\ v\ P, W, Fe, 1;™", T, s} 
the evolution of the fluid elements (Sec 
primitive electron fraction Y^, the conserved electron 
fraction Y^°'^, the temperature T, and the specific 



describing 
(Sec 



IIBl 



of variables 
describing 
The 



entropy s are only necessary when microphysical finite- 
temperature equations of state are used. In addition to 
these two sets of variables, we also need to consider the 
"pseudo-evolved" atmosphere mask Ma described in Ap- 
pendixjC] Thus, in total, we have 24-1-19-1-1 = 44 evolved 
components that potentially need to be synchronized. 

As described in Sec. |IIB| the update terms for the 
spacetime variables are computed via finite differences 
and thus require ghost-zone synchronization after each 
intermediate step. In addition, they are also subject to 
AMR buffer-zone synchronization via prolongation to ob- 
tain valid ghost data from the coarse grid in the buffer 
zone. 



As described in Sec. |II A| the update terms for the 
evolved conserved fluid variables are computed from re- 
constructed primitive variables at cell interfaces and thus 
also require ghost and buffer-zone synchronization in the 
same way as the spacetime variables. The conservative 
to primitive conversion requires the conserved variables 
and valid initial guesses for the primitive variables. Typ- 
ically, these initial guesses are taken from the last valid 
time step on the given cell. Since cells located in the 
buffer zone become invalid during time integration sub- 
steps and need to be refilled via buffer-zone prolongation 
after a full time step, we also need to synchronize those 
primitive variables that are used as initial guesses in the 
conservative to primitive conversion. In our case, these 
are p, e, w*, and T. Note that we do not need to syn- 
chronize the global primitive velocity since it is later 
obtained from a coordinate transformation. 

Furthermore, we need to update the atmosphere mask 
Ma in each intermediate step via inter-processor and 
inter-patch synchronization, and also via buffer-zone pro- 
longation after each full time step. This is necessary be- 
cause the atmosphere mask is only set on cells of the 
evolved grid (i.e. all cells excluding ghost zones). Op- 
erations like conservative to primitive conversion, which 
depend on the atmosphere mask, are performed on the 
entire grid, including ghost zones. Thus they require a 
synchronized atmosphere mask. In addition, the syn- 
chronization order of the atmosphere mask is important 
during buffer-zone prolongation. Before prolongating all 
other required quantities, we first prolongate the atmo- 
sphere mask. Immediately afterwards, cells are set to at- 
mosphere values according to the atmosphere mask. The 
atmosphere mask itself is cleared (also see Appendix [C| . 
This completes the evolution step, and all variables are 
in their final state for the given evolution step. Now, it 
is possible to prolongate also all remaining variables as 
discussed above. 

Finally, we need to synchronize all variables (except 
for the atmosphere mask^) via prolongation after regrid- 



IIAl 



The atmosphere mask does not need to be synchronized because 
it is not valid during regridding. As explained in Appendix[C] it 
is only valid during time integration substeps where regridding 
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TABLE VII: Required synchronizations for each quantity for the three synchronization operations. See text for more details. 



Operation 


inter-processor / inter-patch sync. 


prolongation (buffer zone) 


prolongation (regridding) 


Quantities 












{D,T,S„Yr'',P,e,v\T} 


{D,T,S,,Yr'',P,e,v\Y,,T,s} 




{Ma} 


{Ma} 



ding. A subsequent conservative to primitive conversion 
ensures that the two conservative and primitive sets of 
hydrodynamical variables are consistent with each other. 
Even though regridding requires all variables to be syn- 
chronized and is thus rather expensive, fortunately, this 
operation usually does occur only infrequently, say every 
64 iterations, when moving the fine grids during binary 
neutron star evolution, and only very infrequently, say 
every couple of thousands of iterations, when adding ad- 
ditional refinement levels during stellar collapse or neu- 
tron star collapse. 

In Table |VIH we explicitly list all quantities that must 
be updated during one of the three possible synchro- 
nization operations. The most frequent operation, inter- 
processor and inter-patch synchronization require the 
least number of variables to be updated. Prolongation 
during regridding, which is the least frequent synchro- 
nization operation, requires the full set of variables (ex- 
cept for the atmosphere mask Ma which is invalid outside 
of a full time integration step) . Also note that the global 
primitive velocity never needs to be synchronized be- 
cause it is obtained from the local primitive velocity 
via a coordinate transformation after each synchroniza- 
tion step. Similarly, the Lorentz factor W and the pres- 
sure P are never synchronized since they are computed in 
the conservative to primitive routine, which is exectued 
after each synchronization operation. 



Appendix E: Volume Integration 

Several quantities in our code require volume integra- 
tion over the entire numerical grid. For instance, the 
total baryonic mass is given by 



imated numerically by 



(fix D{x, y, z) 



(El) 



in terms of the conserved density D in the Cartesian ten- 
sor basis*. In Cartesian coordinates, this can be approx- 



is not allowed. We clear it in any new grid region. 
We remark that our code uses the conserved density D in the local 
coordinate basis. Since D is a densitized scalar, | |E1| requires an 
additional Jacobian factor to transform D to the global basis. 
For simplicity of discussion, we omit this here and temporarily 
assume that D is given in the global basis. 



ijk 



(E2) 



where Ax, Ay, and Az is the grid spacing and the indices 
i,j,k, in this context, denote grid indices. In generic 
curvi-linear coordinates, the global grid spacing is not 
constant anymore. In order to compute the volume inte- 
gral with respect to global coordinates, we make use of 
the local volume element 



d M = AuAvAw . 



(E3) 



where Au, Av, and Aw denotes the local uniform grid 
spacing, and we make use of the relation between local 
volume form d^u aird global volume form 



det 



dx^ 



(E4) 



The volume form d^x is introduced as an additioiral grid 
functioir which can be computed once the coordinates 
and grids are set up. 

Next, we need to take into account the non-trivial over- 
lap between neighboring grid patches. For instance, the 
spherical boundary of the spherical outer grid (Fig. [ij 
cuts through cells of the central Cartesian patch, i.e., 
parts of the Cartesian cells reach into the nominal do- 
main of the spherical grid. Consequently, the volume 
associated with each of those cells is only a fraction of 
the volume of the eirtire cell. In practice, we set up a 
weight mask Wijk defining the contribution of each cell 
to the total volume. A cell fully contained on the nom- 
inal grid has a weight of Wijk = 1- Correspondingly, a 
cell completely outside of the nominal grid has a weight 
of Wijk = 0. Cells, whose vertices are not all on the nom- 
inal grid, carry a weight < Wijk < 1- In that case, we 
determine the weight by using 3D Monte-Carlo integra- 
tion [e.g., Illlj of the volume fraction of the overlapping 
regions. The weights need to be calculated only once af- 
ter the grids have beeir setup and therefore the cost of 
Monte Carlo volume integration is negligible compared 
to the total cost of the simulation. 

For simplicity, we absorb the weight mask into the vol- 
ume form (E4), i.e., we effectively store 



{d^x)ijk = AuAvAv 



det 



dx 



ijk ■ 



(E5) 



ijk 
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FIG. 26: Binary neutron stars: influence of the outer bound- 
ary on the accuracy of the wave extraction and evolution. The 
upper panel shows the "+" polarization of the Weyl scalar 
D^4, extracted via CCE for the two setups with different 
outer boundary locations. At time t ~ 7.5 ms, when the 
outer boundary in the setup with Rb ~ 2OOOM0 comes in 
causal contact with the interior evolution, differences start to 
become visible for the Rb ~ 2000Mq setup: the amplitude 
of _D*I'4 deviates by ~ 7%, the phase (j> deviates by ~ 0.2 rad, 
and the Li-norm of the Hamiltonian constraint is larger 

by ~ 15%. 



where the indices i,j,k label grid points and are not sub- 
ject to the Einstein sum convention. Similar to the Jaco- 
bians introduced for computing global Cartesian deriva- 
tives from local finite differences, any volume integration 



needs to take into account ( E5 1 
the form 



For instance ( E2 1 takes 



Mb =^D,jk{d^x)^jk- 

ijk 



(E6) 



Appendix F: Influence of the Outer Boundary 

All GR binary neutron star merger simulations to date 
employ grids which are too small to allow for causally dis- 



connected outer boundaries. Since no constraint preserv- 
ing boundary conditions are known for the BSSN evolu- 
tion system, the simulations may be affected by incoming 
constraint violations. Thus, it is interesting to investigate 
the influence of the outer boundary condition on the inte- 
rior evolution and extracted GWs of the binary neutron 
star merger problem considered in Sec. |IIID| when the 
boundary is not causally disconnected. 

We compare a simulation with outer boundary at 
Rb = 2000AfQ to the simulations in Sec. HID which 
use an outer boundary at Rb — 28OOM0. The setup 
with _Rb = 2OOOM0 has an outer boundary wich is 
in causal contact with the interior evolution and the 
wave-extraction region during the merger and ring-down 
phases. All simulations impose an approximate and 
non-constraint preserving radiative boundary condition 
(e.g. [21]). We focus on baseline resolution rl. We ex- 
pect the simulations to be very similar at least up to 
the point when the constraint violations from the outer 
boundary reach the wave-extraction region which hap- 
pens at i ~ 7.5 ms. 

In Fig. [26] we show the "+" polarization of the leading 
order harmonic (£, m) = (2, 2) mode of the complex Weyl 
scalar 1)^4 computed via CCE. The difference in ampli- 
tude are on the order of ^ 7%. The effects on the phase 
are more subtle and not clearly visible from a simple in- 
spection of the waveform itself. Therefore, in the two 
panels below, we plot the phase of the (£, m) = (2,2) 
mode. The maximum dephasing in the two simulations is 
^ 0.2 rad and thus, the systematic dephasing due to the 
influence from the outer boundary is only slightly below 
the one due to the convergent numerical error. This in- 
dicates that when the resolution is further increased, the 
error due to constraint violations from the outer bound- 
ary cannot be neglected anymore. 

In the same, flgure, we also show the Li-norm^ of the 
Hamiltonian constraint ||i?||i for the two simulations. 
We flnd that the difference of ~ 15% is smaller than the 
difference of ~ 25% between the numerical resolutions rl 
and r2, but not so small that it can be ignored. 

Finally, we also compare mass and spin of the merger 
remnant, and flnd that the differences are on the order 
of the numerical error between resolutions rl and r2. 

Overall, we flnd that causally disconnected outer 
boundaries have a non-negligible impact on the accu- 
racy of the binary neutron star simulation presented in 
Sec. [hid] It is thus likely that longer inspiral simulations 
are even more strongly affected. 
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